arXiv:1506.02466v3 [nucl-th] 27 Jan 2016 


Uncertainty analysis and order-by-order optimization of chiral nuclear interactions 


B. D. Carlsson,^’* A. Ekstrom,^’^’ 't' C. Forssen,^’ ^ D. Fahlin Stromberg,^ 

G. R. Jansen,^’^ O. Lilja/ M. Lindby/ B. A. Mattsson/ and K. A. Wendt^’^ 

^Department of Fundamental Physies, Chalmers University of Teehnology, SE-412 96 Goteborg, Sweden 
^Department of Physies and Astronomy, University of Tennessee, Knoxville, TN 37996, USA 
^Physies Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA 
National Center for Computational Scienees, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA 

(Dated: January 28, 2016) 

Chiral effective field theory (yEFT) provides a systematic approach to describe low-energy nuclear 
forces. Moreover, yEFT is able to provide well-founded estimates of statistical and systematic 
uncertainties — although this unique advantage has not yet been fully exploited. We fill this 
gap by performing an optimization and statistical analysis of all the low-energy constants (LECs) 
up to next-to-next-to-leading order. Our optimization protocol corresponds to a simultaneous fit 
to scattering and bound-state observables in the pion-nucleon, nucleon-nucleon, and few-nucleon 
sectors, thereby utilizing the full model capabilities of yEFT. Finally, we study the effect on other 
observables by demonstrating forward-error-propagation methods that can easily be adopted by 
future works. We employ mathematical optimization and implement automatic differentiation to 
attain efficient and machine-precise first- and second-order derivatives of the objective function with 
respect to the LECs. This is also vital for the regression analysis. We use power-counting arguments 
to estimate the systematic uncertainty that is inherent to yEFT and we construct chiral interactions 
at different orders with quantified uncertainties. Statistical error propagation is compared with 
Monte Carlo sampling showing that statistical errors are in general small compared to systematic 
ones. In conclusion, we find that a simultaneous fit to different sets of data is critical to (i) identify 
the optimal set of LECs, (ii) capture all relevant correlations, (iii) reduce the statistical uncertainty, 
and (iv) attain order-by-order convergence in yEFT. Furthermore, certain systematic uncertainties 
in the few-nucleon sector are shown to get substantially magnified in the many-body sector; in 
particlar when varying the cutoff in the chiral potentials. The methodology and results presented 
in this Paper open a new frontier for uncertainty quantification in ab initio nuclear theory. 


I. INTRODUCTION 

Uncertainty quantification is essential for generating 
new knowledge in scientific studies. This insight is res¬ 
onating also in theoretical disciplines, and forward er¬ 
ror propagation is gaining well-deserved recognition. For 
instance, theoretical error bars have been estimated in 
various fields such as neurodynamics [1], global climate 
models [2], molecular dynamics [3], density functional 
theory [4], and high-energy physics [5]. 

In this paper, we present a systematic and practical ap¬ 
proach for uncertainty quantification in microscopic nu¬ 
clear theory. For the first time, we provide a common 
statistical regression analysis of two key frameworks in 
theoretical nuclear physics: ab initio many-body meth¬ 
ods and chiral effective field theory (yEFT). We supply a 
set of mathematically optimized interaction models with 
known statistical properties so that our results can be 
readily applied by others to explore uncertainties in re¬ 
lated efforts. 

The ab initio methods for solving the many-nucleon 
Schrodinger equation, such as the no-core shell model 
(NCSM) [6] and the coupled cluster (CC) approach [7], 
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are characterized by the use of controlled approxima¬ 
tions. This provides a handle on the error that is 
associated with the solution method itself. Over the 
past decade there has been significant progress in first- 
principles calculations of bound, resonant, and scatter¬ 
ing states in light nuclei [6, 8-12] and medium-mass nu¬ 
clei [7, 13-16]. The appearance of independently con¬ 
firmed and numerically exact solutions to the nuclear 
many-body problem has brought forward the need for an 
optimized nuclear interaction model with high accuracy, 
quantified uncertainties, and predictive capabilities. 

yEFT is a powerful and viable approach for describ¬ 
ing the low-energy interactions between constituent nu¬ 
cleons [17, 18] — a cornerstone for the microscopically 
grounded description of the atomic nucleus and its prop¬ 
erties. Most importantly, the inherent uncertainty of the 
yEFT model can be estimated from the remainder term 
of the underlying momentum-expansion of the effective 
Lagrangian. We refer to this error as a systematic model 
uncertainty. 

We use the common term low-energy constants (LECs) 
to denote the effective parameters of a nuclear interac¬ 
tion model. Indeed, for the description of atomic nuclei, 
the numerical values of the LECs play a decisive role. In 
the yEFT approach, the LECs can in principle be con¬ 
nected to predictions from the underlying theory of quan¬ 
tum chromodynamics (QCD), see e.g. Ref. [19]. How¬ 
ever, the currently viable approach to accurately describe 
atomic nuclei in yEFT requires that the LECs are con- 
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strained from experimental low-energy data. The bulk of 
this fit data traditionally consists of cross sections mea¬ 
sured in nucleon-nucleon (A/TV) scattering experiments. 
Most often, this data is parameterized in terms of phase 
shifts [20, 21]. However, experimental data comes with 
error bars, which implies that a thorough statistical error 
analysis of the constructed nuclear Hamiltonian can only 
be performed when fitting directly to nuclear scattering 
observables. This optimization procedure gives rise to 
statistical uncertainties on the LECs. 

In general, the determination of the LECs constitutes 
an extensive nonlinear optimization problem. That is, 
the relatively large number of parameters makes it chal¬ 
lenging to find optimal values such that the experimental 
fit data is best reproduced. Various methods and objec¬ 
tive functions have been used to solve this problem for 
a wide array of available nuclear-interaction models [22- 
28]. More often than not, the parameters of the mod¬ 
els were fitted by hand. Mathematical optimization al¬ 
gorithms were only recently introduced in this venture 
by Ekstrom et al. [29] and by Navarro Perez et al. [30]. 
Eirst attempts to investigate the statistical constraints on 
the LECs of mathematically optimized interactions have 
recently been performed in the NN sector with coarse¬ 
grained ^-shell interactions [31-34] and with yEET NN 
interactions [35, 36]. 

The so called power-counting scheme of the yEET ap¬ 
proach offers a systematically improvable description of 
TVTV, three-nucleon (TVTVTV), and pion-nucleon (ttTV) in¬ 
teractions. It provides a consistent framework in which 
LECs from the effective ttTV Lagrangian also govern the 
strength of pion-exchanges in the NN potential and of 
long- and intermediate-range TVTVTV forces. This implies 
that ttTV scattering data can be used to constrain some 
LECs that enter the chiral nuclear Hamiltonian. 

Eurthermore, yEET offers an explanation for the ap¬ 
pearance of many-nucleon interactions, such as TVTVTV- 
diagrams, and the fact that they provide higher-order 
corrections in the hierarchy of nuclear forces. Still, effec¬ 
tive TVTVTV forces are known to play a prominent role in 
nuclear physics [8, 37, 38]. Most often, the LECs that are 
associated with the TVTVTV terms have been determined 
relative to existing TVTV Hamiltonians. These LECs are 
optimized against a few select binding energies, excita¬ 
tion energies, or other properties of light nuclei. 

The extended approach that is presented here is con¬ 
ceptually consistent with yEET in the sense that the 
TVTV+TVTVTV Hamiltonian is constrained from a simultane¬ 
ous mathematical optimization to TVTV and ttTV scattering 
data, plus observables from TVTVTV bound states including 
the electroweak process responsible for the /3-decay of ^H. 
Eurthermore, we include the truncation error of the chiral 
expansion to take systematic theoretical errors into ac¬ 
count. If correctly implemented, the truncation error of 
an observable calculated in this scheme should decrease 
systematically with increasing order in the yEET expan¬ 
sion. Indeed, we will show that the resulting propagated 
uncertainties of a simultaneous fit are smaller and exhibit 


a more obvious convergence pattern compared to the tra¬ 
ditional separate or sequential approaches that have been 
published so far. 

Below, we summarize the work presented in this Paper 
by listing three specific objectives: 

• Establish a systematic framework for performing 
mathematical optimization and uncertainty quan¬ 
tification of nuclear forces in the scheme of yEET. 
Our approach relies on the simultaneous optimiza¬ 
tion of the effective nuclear Hamiltonian to low- 
energy ttTV, TVTV, and TVTVTV data with the inclusion 
of experimental as well as theoretical error bars. 

• Demonstrate methods to propagate the statisti¬ 
cal errors in the order-by-order optimized nuclear 
Hamiltonian to various nuclear observables and in¬ 
vestigate the convergence of the chiral expansion. 

• Deliver optimized chiral interactions with well- 
defined uncertainties and thoroughly introduce the 
accompanying methodological development such 
that our results can be easily applied in other cal¬ 
culations. 

Our paper is organized as follows: In Section H we 
introduce the methodology. We start with the construc¬ 
tion of the nuclear potential from yEET and proceed to 
the calculation of observables and the optimization of pa¬ 
rameters. In particular, we introduce automatic differen¬ 
tiation for numerically exact computation of derivatives, 
and we discuss the error budget and error propagation. 
The results of our analysis, for potentials at different 
orders in the chiral expansion and using different opti¬ 
mization strategies, is presented in Section HI. We study 
the order-by-order convergence, the correlation between 
parameters, and we present first results for few-nucleon 
observables with well-quantified statistical errors propa¬ 
gated via chiral interactions. The consequences of our 
findings in the few- and many-body sectors are discussed 
in Section IV, and in Section V we present an outlook for 
further work. 


II. METHOD 

In this section we give an overview of the nuclear 
yEET that we employ to construct a nuclear poten¬ 
tial (Sec. HA). The optimal values for LECs are not 
provided by yEET itself; they need to be constrained 
from a fit to data. Eor completeness we will summa¬ 
rize the well-known methods to calculate the relevant 
experimental observables: TVTV scattering cross sections 
(Sec. HB), TVTV effective range parameters (Sec. HC) 
and bound state properties for A < 4 nuclei using the 
Jacobi-coordinate no-core shell model (Sec. HD). We also 
present the objective function (Sec. HE), the optimiza¬ 
tion algorithm (Sec. HE), and the formalism for the sta¬ 
tistical regression (Sec. HG). 
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A. The nuclear potential from xEFT 


NN NNN 


The long-range part of the nuclear interaction in xEFT 
is governed by the spontaneously broken chiral symmetry 
of QCD and mediated by the corresponding Goldstone- 
boson; the pion (tt). This groundbreaking insight [39] en¬ 
ables a perturbative approach to the description of phe¬ 
nomena in low-energy nuclear physics [40]. High-energy 
physics that is not explicitly important is accounted for 
through a process of renormalization and regularization 
with an accompanying power counting scheme. The 
expansion parameter is defined as Q/A^, where Q is 
associated with the external momenta (soft scale) and 
Mp (hard scale), with Mp ^ 800MeV the mass 
of the rho meson. The benefit of a small-parameter ex¬ 
pansion is that higher orders contribute less than lower 
orders. If the series is converging, an estimate of the 
magnitude of the truncation error is given by the size of 
the remainder. 

The chiral order of a Feynman diagram is governed by 
the adopted power-counting scheme. Given this, any chi¬ 
ral order z/ > 0 in the expansion will be identified with a 
finite set of terms proportional to {Q/. In this work 
we have adopted the standard Weinberg power-counting 
(WPG) which is obtained from the assumptions of naive 
dimensional analysis. For the scattering of two or more 
nucleons without spectator particles, v is determined by 
(see e.g. Ref. [18]) 



NNLO 

u = 3 




Figure 1. Schematic overview of the Feynman diagrams 
present at leading order (LO), next-to-leading order (NLO), 
and next-to-next-to-leading order (NNLO). Nucleons (pions) 
are represented by solid (dashed) lines. The three-nucleon 
(NNN) interaction enters at NNLO. A circle, diamond and 
square represents a vertex of order A = 0, 1 and 2 respec¬ 
tively. 


u = 2A-4 + 2L + '^Ai (1) 

where A is the number of nucleons and L is the number 
of pion loops involved. The sum runs over all vertices 
i of the considered diagrams and is proportional to 
the number of nucleon fields and pion-mass derivatives 
of vertex i. A^ > 0 for all diagrams allowed by chiral 
symmetry. In Fig. 1 we show the different interaction di¬ 
agrams that enter at various orders. For the NN system, 
contributions at = 1 vanish due to parity and time- 
reversal invariance. Also, we consider nucleons and pions 
as the only effective degrees of freedom and ignore pos¬ 
sible nucleon excitations, i.e., we use the so-called delta¬ 
less version of yFFT. 

The interaction due to short-range physics is parame¬ 
terized by contact terms, which also serve to renormalize 
the infinities of the pion loop integrals. The order-by¬ 
order expansion of this zero-range contribution is also 
organized in terms of increasing powers of Q/A^. Due to 
parity, only even powers of u are non-zero. Furthermore, 
the contact terms of order z^ = 0 contribute only to par¬ 
tial waves with angular momentum L = 0, i.e. S'—waves, 
whereas v = 2 contact terms contribute up to P-waves. 
In general, the contact interaction at order v acts in par¬ 
tial waves with L < vj2. Following Fq. (1), the terms 
in the yEFT expansion, up to third order, are given by 
a sum of contact interactions Wt and one- plus two-pion 


exchanges, denoted by and V 27 r 5 respectively: 

uo = vf + 

w = vlo++ k'^ (2) 

VnNLO = UfLO + vk + kw^ + VnNN- 

The superscript indicates the separate chiral orders u = 
0,2,3, referred to as leading-order (LO), next-to-leading 
order (NLO), and next-to-next-to-leading order (NNLO). 
For detailed expressions, see e.g. Ref. [27]. The three- 
nucleon interaction, Vnnn^ contains three different dia¬ 
grams as shown in Fig. 1. These correspond to two-pion 
exchange, one-pion exchange plus contact, and a pure 
NNN contact term. Insofar, the analytical expressions 
for the NN potential have been derived up to fifth or¬ 
der (N4LO) [41, 42]. The partial-wave decomposition for 
the NNN interaction at NNLO is well known [43], while 
the N3LO contribution was published very recently [44]. 
Note that the connected four-nucleon diagrams also ap¬ 
pear at this higher order. In the present work we limit 
ourselves to NNLO for completeness. 

The strengths of the terms in the yEFT interaction 
are governed by a set of LECs. These parameters play 
a central role in this work, and we discuss in detail how 
they are constrained from measured data. In general, for 
each chiral order there will appear a new set of LECs. For 
the nuclear interactions used in this work, see Eq. (2), the 
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corresponding LECs are denoted 

, Cb, , C3p„ , C.p ,, Cap,, CapJ, 

{Cl,C3,C4}, 

VnNN ^ {ci,Cs,C4^,Cd,Ce}- 


Furthermore, there are additional constants that must 
be determined before making quantitative predictions in 
XEFT. Here, we set the axial-vector coupling constant to 
the experimentally determined value of qa = 1.276 [45] 
for LO, whereas for the higher orders we use the renor¬ 
malized value of gA = 1-29 to account for the Goldberger- 
Treiman discrepancy [27]. At all orders we use = 
92.4 MeV [27]. All other physical constants, such as nu¬ 
cleon masses and the electric charge, are taken from CO¬ 
DATA 2010 [46], except the pion masses for which we 
have used the values from the Particle Data Group [47]. 

Note that LECs that determine the sub-leading irN in¬ 
teraction vertices occur in both the NN interaction and 
the two-pion-exchange part of the ATVTV, see Refs. [27, 
43]. Besides offering this pion-vertex link between the 
NN and the NNN interaction, the irN interaction model 
of yEFT allows to describe irN scattering processes. 
Consequently, experimental ttN scattering data can be 
used to constrain the long-range part of the nuclear in¬ 
teraction. The lowest order terms of the effective irN 
Lagrangian have u = 1 and are free from LECs, besides 
gA and F^. At order v = 2 the LECs ci, C 2 , C 3 , and C 4 
enter. Higher-order eN LECs, such as di + (i 2 , ds, ds 
and di 4 — di 5 , enter at z/ = 3 while 614 to eig appear at 
z/ = 4. In total, there are 13 LECs in the ttN Lagrangian 
up to fourth order. 

The different masses and charges of the up and down 
quarks give rise to isospin-violating effects [18, 27]. There 
are both short- and long-range isospin-violating effects. 
The long-range effects are of electromagnetic (EM) origin 
and for this contribution we use the well-known set of 
potentials 


=^C1 + ^C2 + Wp + k[M^ 




(np) 


EM 


x.(np) 

■’^MM ’ 


( 4 ) 


where Cl is the static Coulomb potential, C 2 the rel¬ 
ativistic correction to the Coulomb potential [48], VP 
is the vacuum polarization potential [49], and MM the 
magnetic-moment interaction [50]. The long-range ef¬ 
fects become increasingly important as the scattering en¬ 
ergy approaches zero; consequently we include all the 
above long-range effects at all orders in the chiral ex¬ 
pansion. We also consider short-range isospin-breaking 
mechanisms. At NLO, the Cisq contact is split into three 

charge-dependent terms: and . At this 

order, and above, we also take the pion-mass splitting 
into account in one-pion exchange terms [18]. 

An effective field theory often has to handle more than 
one expansion parameter. In our case, the nucleon mass. 


Mat = 2MpMn/{Mp + M^) where Mp (M^) is the pro¬ 
ton (neutron) mass, provides such an extra scale and the 
use of the heavy-baryon chiral perturbation theory intro¬ 
duces relativistic corrections with factors of I/Mat- We 
count these corrections as Q/M^ ~ (Q/A^)^ [40, 51]. 
This choice implies that no relativistic corrections ap¬ 
pear in the NN sector up to the order considered in this 
paper. 

To regularize the loop integrals that are present in the 
two-pion exchange diagrams we employ spectral func¬ 
tion regularization (SFR) [52] with an energy cutoff 
A = 700 MeV. The nuclear interaction is calculated per- 
turbatively in yEFT. A nuclear potential that can be 
used for bound and scattering states is obtained by iter¬ 
ating the terms of the chiral expansion in the Lippmann- 
Schwinger or Schrodinger equation [53]. We employ the 
minimal-relativity prescription from Ref. [54] to obtain 
relativistically-invariant potential amplitudes. The ul¬ 
traviolet divergent Lippmann-Schwinger equation also re¬ 
quire regularization. We remove high-momentum contri¬ 
butions beyond a cutoff energy A by multiplying the NN 
and NNN interaction terms with standard (non-local) 
regulator functions fNN{p) and fNNN{p^Q)^ respectively. 


/jvjv(p) = exp 



( 5 ) 


and 


fNNN{p, q) = exp 


/ y 

V 4A2 ) J ’ 


(6) 


where p and q are the Jacobi momenta of the interact¬ 
ing nucleons. In this work, we mainly use A = 500 MeV 
and n = 3. However, we also explore the consequences 
of varying A in steps of 25 MeV between 450 — 600 MeV. 
The canonical power-counting, i.e. WPC, and the non- 
perturbative renormalization of nuclear yEFT in its cur¬ 
rent inception is currently under some debate [55, 56]. In 
relation to this it should be stressed that our implementa¬ 
tion of statistical regression methods and gradient-based 
optimization methods furnishes an independent frame¬ 
work to extract well-founded estimates of the uncertain¬ 
ties in theoretical few-nucleon physics and a tool to assess 
the convergence properties of yEFT. 


B. Nuclear scattering 

The NN scattering observables are calculated from the 
spin-scattering matrix M [57, 58]. This is a 4 x 4 matrix 
in spin-space that operates on the initial state to give the 
scattered part of the final state. Thus, M is related to 
the conventional scattering matrix S' by M = |^(S — 1), 
where p is the relative momentum between the nucleons. 
The decomposition of M into partial waves is given by 
(see e.g. [ 20 ]) 
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M* f 

±v±m'm 






L' 


J 


s J 
m —m 


{L\s'\S^ -l\L,s), (7) 


where the big parentheses are Wigner 3j-symbols, s {s') and m {m') are initial (final) total spin and spin projection, 
respectively, J = L + s is the total relative angular momentum and L (J) is 21/ + 1 (2 J + 1). The quantization axis is 
taken along the direction of the incoming nucleon and 0 gives the center-of-mass scattering angle. The /S-matrix for 
the scattering channel with angular momentum J can be parameterized by the Stapp phase shifts [59], 

„j _f cos2ej sin2ej\ . . 

L=j±i - cos2ej ) 


for the coupled triplet channel, and 


S' 


J _ 
L=J — 


/ ^2i6j 2yj 

sin 27 j 


^g*(^j+^j,j) 27 j\ 

cos 27 j j 


(9) 


for the (coupled) singlet-triplet channel with L = J. 


The spin-singlet {S = 0) phase shift is denoted by 6l=j^ 
the spin-triplet {S = 1 ) phase shift by Sl^j^ while ej 
represents the triplet-channel mixing angle and 7 j is the 
spin-flip mixing angle [60] ( 7 j = 0 for pp scattering). 

In practice, the infinite sums in Eq. (7) are truncated 
at I/,I/' < I/max- Calculations that involve long-ranged 
EM effects require Tmax > 1000 in order to reach conver¬ 
gence, while I/max = 30 is sufficient for the part coming 
from the short-ranged nuclear interaction. This leads 
to a natural separation of the terms in Eq. (7), see e.g. 
Ref. [50]. In brief, all EM amplitudes are calculated in¬ 
dependently in Coulomb Distorted-Wave Born Approx¬ 
imation (CDWBA) using Vincent-Phatak matching [61] 
to handle the difficulties of the Coulomb interaction in 
momentum space. Eor the channel, the C 2 and VP 
interactions are strong enough that a small correction to 
the bare phase shifts is needed, resulting in 


^total — ^ 


(CDWBA) 

Cl+ATAT 


+ Ao - po - To 


( 10 ) 


where phase shift of the Coulomb and 

the chiral NN interactions computed in CDWBA, po (tq) 
is the C 2 (VP) phase shifts in CDWBA, and Aq is a 
correction calculated by interpolating between the values 
tabulated by Bergervoet et al [62]. In principle, Aq is 
dependent on the interaction model for the strong force; 
this effect has been shown to be very small [62] and was 
not considered here. 

We compute the VP phase shifts, r^, in CDWBA us¬ 
ing the variable-phase method [63] . The values we obtain 
agree with the ones that are tabulated by Bergervoet 
et al. [62]. The VP amplitude is calculated in the first- 
order approximation derived by Durand [49] using the ex¬ 
pansion parameter X = {T\g,}jMp{l — cos{0)))^ where 

rUe is the electron mass. We find that X < 0.031 for all 
scattering data that is employed in this work. The MM 


amplitude for np and pp scattering is given by Stoks^ [50] . 

The Stapp phase shifts are calculated from the real¬ 
valued free reaction matrix R [64], which is defined 
through a Lippman-Schwinger type equation [64] 


L^P ^P) ~ ^L'L (P ’ P) 


r,, Jo 


P P _ rrP ’ 




( 11 ) 


where V is the potential, p the reduced mass, and V 
denotes the Cauchy principal value. 

Due to parity and time-reversal invariance, the scatter¬ 
ing matrix M has six linearly independent elements. We 
employ the Saclay parameterization [57], with complex 
amplitudes a to /, to express 

M(q, k) = i {(a + 5 ) -1- (a - b)ai • ra2 • r 

+ (c + d){(Ti ■ q)(cr2 ■ q) + (c - d){(Ti ■ k){(T2 ■ k) 

- e{(Ti + 0-2) ■ f - f{(Ti - 0-2) ■ r}, 

where q = p' — p is the momentum transfer, k = (p' + 
p )/2 and r = q x k. Eor identical particles, / will be 
zero. Expressions for the scattering observables in terms 
of the Saclay parameters can be found in Ref. [57] for 
identical particles and in Ref. [58] for the more general 
case of non-identical particles. 

Eor the theoretical description of the irN scattering 
observables we use the fourth order yEET expressions 
according to Refs. [65, 66]. A detailed description of the 
EM amplitudes that we employ are given in Refs. [67-70] 


^ Note that Eq. (24) in Ref. [50] has the wrong sign. Furthermore, 
Eq. (25) should have |sin(^)|. 
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C. Effective range parameters 

The effective-range expansion (ERE) of low-energy 
phase shifts [71] provides parameters that can be directly 
compared to experimentally inferred values. The ERE 
can be expressed in the general form 

A{p) + = - ^ + y V + 0(/). 

(13) 

The functions A{p) and B{p) depend on the choice of 
included long-range EM effects and ^^lr+atat is the phase 
shift of the total nuclear potential (long-range plus strong 
NN) relative to the phase shift of only the long-range 
part. 

Eor nn and np scattering we have A{p) = 0 and B{p) = 
1 [71] since there are no EM effects. The corresponding 
ERE parameters are denoted a^p and r^p. 

Eor pp scattering we calculate ERE parameters for 
the nuclear plus Coulomb potential, i.e., using the phase 
shifts expressions for Ac{p) and Bc{p) 

can be found in Refs. [62, 71]. The corresponding ERE 
parameters are denoted and r^. 

In practice, the ERE parameters are determined using 
a linear least-squares fit to 20 equally-spaced phase shifts 
in the Tiab = 10 — 100 keV range. 


D. Few-nucleon observables 

We employ the Jacobi-coordinate version of the 
NCSM [72] to compute bound-state observables for 
and ^’^He. Apart from binding energies and radii we also 
compute the deuteron quadrupole moment, (5(^H), and 
the comparative half-life for the triton /Ti/ 2 (^H). 

In the NCSM, observables and wave functions are ob¬ 
tained from the exact solution of the eigenvalue problem 
B I'lp) = E . In this work, the nuclear Hamiltonian B 
is given by 


AAA 

^ = E + E R + E (14) 

i<j=l i<j=l i<j<k=l 


where Tij are relative kinetic energies while Vij and Vijk 
are the NN and NNN interactions, respectively. In our 
calculations we use the isoscalar approximation as pre¬ 
sented in Ref. [73]. The model-space dimension is deter¬ 
mined from the maximal number of allowed harmonic- 
oscillator (HO) excitations A^max- We obtain essentially 
converged results in a HO basis with oscillator energy 
hw = 36 MeV and model-space dimension A^max = 40(20) 
for A = 3(4). 

The experimentally measured electric-charge radius 
can be related to the theoretically calculated point- 
proton radius through the relation [74] 


2 2 2 ^22 \ 2 
^pt-p ^ch ^DF •) 


(15) 


where (r^) is the proton (neutron) charge mean- 
squared radius and E {N) is proton (neutron) num¬ 
ber. Eurthermore, r^p = is the Darwin-Eoldy cor¬ 
rection [75] and Ar^ includes effects of two-body cur¬ 
rents and further relativistic corrections. We use Vp = 
0.8783(86) fm and = —0.1149(27) fm^ [76]. Eor all nu¬ 
clei, we use Ar^ = 0. 

Precise results for electroweak observables depend on 
two-body nuclear currents and relativistic effects. yEET 
provides a consistent framework for including such cor¬ 
rections and for deriving quantum-mechanical currents, 
such as the electroweak one, from the same Lagrangian as 
the nuclear force. We follow the approach by Gazit et al. 
[77] and compute the triton half-life from the reduced 
matrix element for the J = 1 electric multipole of 
the axial-vector current 

{E^)^\{^Re\\E^fR)\. (16) 

This matrix element is proportional to cd , the LEG that 
also determines the strength of the NN — irN diagram 
of the NNN interaction. As a consequence, the tri¬ 
ton half-life provides a further constraint of the nuclear 
force. The experimentally determined comparative half- 
life, fTij 2 = 1129.6±3s [78], leads to an empirical value 
for (E^) = 0.6848 ± 0.0011 [77]. 

Eor the deuteron quadrupole moment we choose, in¬ 
stead, to fit to the theoretical value obtained from 
the high-precision meson-exchange NN model GD-Bonn, 
Qd = 0-27 [24], with a 4% error bar that more than 
well covers the spread in values using other NN poten¬ 
tial models [18]. 


E. Objective function 


Using the methods to compute observables outlined 
above, the vector a of numerical values for the LEGs at 
a given order in yEET is constrained using experimental 
data. This is accomplished by minimizing an objective 
function defined as 


X' 




ieM 


(a) - 0“py _ 


) =E^iA), (17) 

^ ieM 


where and denote the theoretical and experi¬ 

mental values of observable Oi in the pool of fit data M, 
and the total uncertainty determines the weight of the 
residual, r^. The optimal set of LEGs is defined from 

OL^ = argminy^(a) (18) 

a 


We wish to explore the physics capabilities and limi¬ 
tations of nuclear yEET by forming different objective 
functions and subsequently probing the precision and ac¬ 
curacy of each one in a statistical regression analysis [79]. 
At each chiral order (LO, NLO, or NNLO) we compare 
two different strategies of minimization: simultaneous 
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(sim) and separate (sep). In the “separate” approach 
we first optimize the sub-leading iiN LECs {ci^di^ e^) us¬ 
ing ttN data. Subsequently, we optimize the NN contact 
potential of the nuclear interaction using NN scatter¬ 
ing data, and finally (at NNLO) the NNN interaction 
is determined by fitting cd and ce to the known bind¬ 
ing energies and radii of and ^He, and the compara¬ 
tive /3-decay half life of Besides the first-ever appli¬ 
cation of novel derivative-based optimization techniques 
to this problem, the “separate” approach is very similar 
to the conventional procedure to constrain the descrip¬ 
tion of the nuclear interaction. In contrast, with the “si¬ 
multaneous” approach we optimize all the LECs up to a 
specific-order in yEET at the same time with respect to 
NN and irN scattering data as well as experimentally de¬ 
termined bound-state observables in the two- and three- 
nucleon systems: and ^He. At LO and NLO, the NN 

interaction does not involve any sub-leading ttN ampli¬ 
tudes, nor are there any NNN force terms. Therefore, at 
these orders the sim-potentials are optimized using only 
NN scattering data and the binding energy, radius, and 
quadrupole moment of the deuteron. A summary of the 
data types that were included in the objective function 
for each potential is given in Table I. 


Table I. Objective functions for the various nuclear inter¬ 
actions in this work. Included data types are marked with 
’X’. For sequential optimization, the subscript ’i’ indicates at 
what stage the model is optimized to that data. Excluded 
data-types are indicated with 


Potential 

Scattering data 
NN ttN 

nn ERE 
parameters 

bound-state data 

LOsep 

X 

- 

- 

- 

- 

LOsim 

X 

- 

- 

X 

- 

NLOsep 

Xi 

- 

X2 

- 

- 

NLOsim 

X 

- 

X 

X 

- 

NNLOsep 

X2 

Xi 

- 

- 

X3 

NNLOsim 

X 

X 

- 

X 

X 


The bulk of the experimental data consists of NN and 
ttN scattering cross sections. Eor the NN data we take 
the SM99 database [21] entries with laboratory scatter¬ 
ing energies < 290 MeV, i.e. the pion-production 

threshold, which constitutes a natural limit of appli¬ 
cability for yEET. This results in N^l = 2045 and 
^data ~ data points, including normalization data. 
The number of normalization constants are A^nSrm = 124 
and A^norm = 148. However, we also explore the con¬ 
sequences of varying between 125-290 MeV. Un¬ 
less otherwise stated, our canonical choice is = 

290 MeV and A = 500 MeV. As there is no neutron- 
neutron scattering data, we use the neutron-neutron ^Sq 
scattering length = —18.95(40) fm [18] and effective 
range = 2.75(11) fm [80] to constrain the parameter 
at order NLO. Eor the ttN scattering observables 
we employ the database from the Washington Institute 
group [81], here referred to as the WI08 database. The 


7tN data consists mainly of differential cross sections and 
some singly-polarized differential cross sections for the 
processes tt^ + p ^ tt^ + p and 7r“ + p ^ tt^ + n. Un¬ 
fortunately, the WI08 database contains very little data 
at low scattering energies, which would have been pre¬ 
ferred to constrain the low-energy theory of yEET. In 
fact, there is no scattering data below Tiab = 10.6 MeV. 
Eor this reason, we include all data up to lab energy 
Tlab = 70 MeV and keep all terms up to, and including, 
u = 4 when calculating ttN observables. A lower chi¬ 
ral order does not give a reasonable description of the 
data. This results in N^l^J = 1347 data points includ¬ 
ing Nno^Ji = 110 normalization data. At the optimum, 
it is usually assumed that the residuals are normally dis¬ 
tributed, and that they are all independent of each other. 
If so, then y^(a^) will comply with a chi-squared distri¬ 
bution with Nm - A^norm “ A^« = Aedf “ A^« = A^dof de¬ 
grees of freedom, where N^ denotes the number of LECs 
(i.e., the number of model parameters). In turn, this 
allows for a standard regression analysis. These rather 
strong assumptions of both the model and the data are 
only approximately fulfilled, mainly due to the inherent 
systematic error in yEET. 

The distribution of residuals, r^, for the NNLOsim po¬ 
tential, which will be thoroughly introduced in Sec. Ill A, 
is shown in Eig. 2. It is clear that the residuals are 
not entirely normally distributed, with a skewness of 
—0.38(3) and excess kurtosis of 5.39(6). The main rea¬ 
son for this deviation can be traced to the inclusion of 
a systematic error in the fit. This can produce a con¬ 
sistent over- or underestimate of observables, resulting 
in a non-zero skewness. A non-zero excess kurtosis indi¬ 
cates that the model error sometimes overestimates the 
uncertainty and in other cases underestimates it, causing 
a too sharp peak near zero in the histogram in Eig. 2. 
We stress that the deviations from normality does not 
invalidate the use of y^(a) as an objective function to fit 
the parameters; it just indicates that the minimizer 
will not be a maximum-likelihood estimator. In fact, we 
find that when optimizing NNLOsim using NN scatter¬ 
ing data up to 125 MeV only, to avoid large model errors, 
the skewness and excess kurtosis of the NN scattering 
residuals are significantly reduced; —0.01(6) and 0.6(1), 
respectively. Still, the propagated uncertainties are very 
similar in these two cases Thus, the minimization and 
subsequent regression analysis of the y^(a) function will 
provide valuable insights into both the model and the 
data [79]. 


1. Total error budget 

Eor each residual, the total uncertainty is divided 
into an experimental part and a theoretical part 


2 2 
^ =^exp 


■ CFi 


theo 


^22 2 2 ( 19 ) 

^exp ^numerical ^method ^model 








Residual 

Figure 2. Residual distribution for the NNLOsim potential, 
with a sample mean and standard deviation of —0.04(1) and 
0.977(9), respectively. The deviations from normality, as dis¬ 
cussed in the text, are mainly due to the model error of xEFT. 


The experimental uncertainty (statistical or systematic) 
is provided by the experimenter. Here, we focus on es¬ 
timating the theoretical uncertainty. As a first step, we 
identify three different components: (1) the numerical 
error originating in finite computational precision, (2) 
the method error due to mathematical approximations in 
the solution of the bound-state or scattering problem, (3) 
the model error that is inherent to the truncation of the 
momentum expansion in xEFT. 

The numerical error is the smallest one and several 
new technical developments, such as automatic differen¬ 
tiation for computing derivatives, allow us to generally 
ignore cTnumericai • Howcvcr, some elements of the statis¬ 
tical analysis can potentially become numerically unsta¬ 
ble if the relative errors are too small. In particular, 
this concerns the computation of the covariance matrix 
through the inversion of the Hessian (33). For this rea¬ 
son we impose a minimum relative uncertainty of 0.01%. 
In practice this requirement only affects the error of the 
deuteron binding energy. 

Regarding the method error, the only significant con¬ 
tributions come from truncating the NCSM model space 
and from the use of the isoscalar approximation in calcu¬ 
lations of bound-state observables. Indeed, for all scat¬ 
tering cross sections we include sufficiently many par¬ 
tial waves to construct an exact scattering matrix. We 
estimate the method error of the NCSM calculations 
using a simple exponential extrapolation, F^(A^max) = 
Eoo + ttexp(—6A^rnax), foi* ^ range of different xEFT po¬ 
tentials. However, the uncertainties from the isoscalar 
approximation dominate the truncation error by an or¬ 
der of magnitude. We therefore use the uncertainties 
presented in Ref. [73] as our method error. 

In practice, we combine the method errors with the 
experimental ones to obtain the resulting weight of each 
bound-state observable in the optimization, see Table IT 
In certain cases, the method error is comparative to, or 
larger than, the experimental error. 


Table II. Experimentally determined values and uncertainties 
for ground-state energies (in MeV) and radii (in fm) for ^’^H 
and ^’^He. The quadrupole moment Q(^H) of the deuteron 
is given in fm^ and Ei denotes the reduced transition ma¬ 
trix element related to the /3-decay of ^H. The last column 
is the combined experimental and method errors. For the 
ground-state energies the method error is much larger than 
the experimental one. Table I indicates which observables are 
included in the optimization. Note that the ^He properties 
are not included in the objective function. 



Exp. value 

Ref. 

e’exp+method 

eCr) 

-2.22456627(46) 

[46] 

0.22 X 10“^ 

eCr) 

-8.4817987(25) 

[46] 

0.028 

E^Re) 

-7.7179898(24) 

[46] 

0.019 

SC^He) 

-28.2956099(11) 

[46] 

0.11 


1.97559(78)^ 

[76, 82] 

0.79 X 10“^ 

^pt-p("H) 

1.587(41) 

[76] 

0.041 

»-pt-p("He) 

1.7659(54) 

[76] 

0.013 

»-pt-p("He) 

1.4552(62) 

[76] 

0.0071 

Q^H) 

0.27(1)*’ 


0.01 

E\eR) 

0.6848(11) 

[77] 

0.0011 


^ The experimental value is r\(^Yi) — r^, we still use the 
value of Tn from Ref. [76] 

^ This is not an empirical value, see the text for details. 


The model errors can be labeled as systematic and are 
the most difficult to assess. We follow the most naive 
xEFT estimate and associate a truncation error with the 
effect of excluded higher-order Feynman diagrams. The 
xEFT expansion up to a given chiral order n includes 
all diagrams that scale as {Q/A^Y where Q G 
The remainder of the diagrams could a priori be assumed 
proportional to (Q/A^)^^^. 

For bound-state properties it is not straightforward to 
associate a relevant and system-dependent momentum 
scale; therefore, we will not include systematic theoreti¬ 
cal errors for these observables. Scattering observables, 
on the other hand, have a well-defined center-of-mass 
momentum. As described in section HB, M-matrix ele¬ 
ments are the fundamental quantities that are needed to 
calculate NN scattering observables and can be parame¬ 
terized by the complex-valued Saclay amplitudes a to f. 
Similarly, the ttN non-spin-flip and spin-flip amplitudes 
and determine the ttN scattering observables [66]. 
Therefore, from the above scaling argument we introduce 
a model error in the scattering amplitudes of the form 

ATdeb = , X e {NN, ttN}, (20) 

where Cnn and are two overall constants that need 
to be determined. 

We assume that both the real and the imaginary parts 
of the Saclay amplitudes a — e scale in this manner. The 
nuclear force does not contribute to the / amplitude so 
we do not impose a model error in that amplitude. Since 
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Figure 3. Saclay amplitudes a to e at Ocm = 45° scattering 
angle for the potential NNLOsim. The model error bands 
were extracted according to the discussion in the text. 


the order of magnitude of each scattering amplitude is 
the same, we assign the same constant of proportionality 
to all of them, see e.g. Fig. 3. The same argument applies 
to the ttN amplitudes. 

We set Q = p to capture the increasing uncertainty 
in the model as the energy increases. The definition 
Q = max{p, [17] seems to have a comparatively 
small impact on the theoretical predictions of the model 
as discussed further in Sec. IIIC. 

To determine Cnn and we use the statistical guid¬ 
ing principle that x^/A^dof for both NN and ttN scatter¬ 
ing should be 1 if the objective function x^ follows a chi- 
squared distribution and all errors have been correctly 
accounted for. This leads to an iterative process where 
first the constants are updated, then the LECs are 
optimized using the previously determined Cx, and so on 
until the values of the constants have stabilized. This 
usually requires no more than three iterations. 


F. Optimization algorithms 

The minimization of Eq. (17), is a non¬ 

linear optimization problem. In this work we have em¬ 
ployed three different non-linear least-squares minimiza¬ 
tion methods at different stages during the optimization: 
POUNDerS [83], Levenberg-Marquardt (LM) and New¬ 
ton’s method. POUNDerS is part of the TAO pack¬ 
age [84] and is a so-called derivative-free method. As 
the label indicates, it does not require the computation 
of any derivatives. This makes it very attractive for use 
with applications where differentiation is a formidable 
task; e.g. nuclear energy density optimization [85] and 
previous optimizations of chiral interactions [29, 35, 86]. 
However, in this work, we have managed to make sig¬ 
nificant progress in the optimization problem by imple¬ 
menting automatic differentiation, which enables us to 
extract machine-precise derivatives of the objective func¬ 
tion. Consequently, the whole class of derivative-based 


optimization algorithms becomes readily available. The 
convergence rate is increased considerably with the LM 
method that employs first-order derivatives of the resid¬ 
uals with respect to the LECs. A further improvement 
can be achieved with Newton’s method that uses also 
the second-order derivatives. At LO, the presence of 
only two LECs to parameterize the potential makes it 
a trivial task to minimize the corresponding objective 
functions. However, already at the next order, NLO, the 
optimization requires quite an effort. There are 11 LECs, 
and in order to provide a reasonable start vector ao of 
numerical values for these we make an initial fit to the 
NN scattering phase-shifts published by the Nijmegen 
group [20]. At NNLO there is a total of 26 LECs, since 
we also need to include all the 13 irN LECs up to order 
z/ = 4. Also at this order we carry out an initial fit to NN 
phase-shifts before proceeding with the optimization of 
the complete objective function. The optimization with 
respect to scattering observables in the irN sector could 
proceed without any fits to phase shifts. 

There is always a risk of getting trapped in local min¬ 
ima and the success of the minimization strongly de¬ 
pends on the starting point ao- Extensive searches were 
performed to search for a global minimum, which is de¬ 
scribed in more detail in Sec. HI A. 


1. Automatic differentiation 


Eirst- and second-order derivatives of x^(<^) with re¬ 
spect to the LECs are needed during the minimization 
process and the subsequent statistical regression analy¬ 
sis, i.e., we need to compute 


g0(theo)(^) 

dcXrn 

dOLm^^n 


\/ i^m 
V i, m, n. 


( 21 ) 


The straightforward numerical approach is to approx¬ 
imate the nth-order derivatives with finite differences. 
The general idea is to form appropriate linear combi¬ 
nations of M function evaluations in the vicinity of the 
point of interest. There are, however, a number of issues 
with this method. Eirst, it is prone to large numerical 
errors since differences of large, almost equal, numbers 
are needed. Second, the result can be very sensitive to 
the choice of step size. Eurthermore, it is also a computa¬ 
tionally demanding method since the number of required 
function evaluations grows quickly with the number of 
dependent variables and order of the derivative. Eor in¬ 
stance, a third-order, finite-difference calculation of first 
and second derivatives with respect to all 26 LECs re¬ 
quires M = 3653 function evaluations. Eor these reasons, 
we abandon finite-difference methods and employ instead 
forward-mode automatic differentiation (AD). The basic 
idea of AD is the following: A computer implementation 
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To solve for the two-nucleon i?-matrix (11) at a given 
on-shell scattering energy we use the well-known method 
of Ref. [88]. It recasts the Lippmann-Schwinger equation 
into a matrix equation 

(/ + VZ)R = V, (23) 

where / is the identity matrix, V is the two-nucleon po¬ 
tential, and Z is a simple diagonal matrix defined in 
Ref. [88] . The R-matrix is easily obtained after inverting 
(I -hVZ) using e.g. LU factorization. First- and second- 
order derivatives of the R-matrix with respect to LECs 
a^c and ay are easily obtained using the AD technology 
and the same LU factorization. 


Figure 4. Comparison between calculated first and second 
derivatives of an objective function using finite differences 
(third order) with different step sizes (filled lines) and auto¬ 
matic differentiation (dashed lines). The calculation is done 
at a minimum where the first derivatives should be approxi¬ 
mately zero. Due to cancellation effects the finite-difference 
method cannot correctly reproduce the low values of the 
derivatives for any step size. 


for calculating the observables, or any computational al¬ 
gorithm for that matter, will consist of a chain of simple 
(or intrinsic) mathematical operations; e.g. addition and 
multiplication, elementary functions such as sin and exp, 
and matrix operations. Therefore, by repeatedly employ¬ 
ing the chain rule^ derivatives with respect to the LECs 
can be calculated alongside the usual function evalua¬ 
tions. Using AD, the derivatives of Eq. (21) can actually 
be computed to machine precision, which is far beyond 
the precision of any reasonable finite-difference scheme. 
This accomplishment is illustrated in Eig. 4, where also 
the dependence on the step size for the finite difference 
method is shown for comparison. 

We implement forward-mode AD using the Rapso- 
dia computational library [87]. Eor the calculation of 
first and second derivatives with respect to Noc different 
LECs, Rapsodia requires a total of 


M = 2 




( 22 ) 


derivative calculations. Eor Noc = 26, this results in 
M = 702, thus considerably more efficient than the finite- 
difference approach. Eurthermore, all calculations that 
do not depend on the LECs are performed only once, 
compared to the brute-force implementation of the finite- 
difference scheme that requires a full calculation for every 
function evaluation. Eurthermore, since all LECs enter 
linearly in the momentum-space formulation of the chi¬ 
ral potential it is very easy to calculate the derivatives of 
the potential with respect to the LECs. Thus, the only 
workhorses in our calculations are the R-matrix evalua¬ 
tion (matrix inversion) of the scattering process and the 
solution to the NCSM eigenvalue problem (matrix diag- 
onalization) as we will discuss next. 




(/ + VZ) 


d^R 


dV 


d'^V 

daxda. 


-(/. 


dV dR 

- Zj 


ZR) 

dV 


OR 


da-r da,, da,, da^ 


(24) 

(25) 


We also use the fact that many derivatives are exactly 
zero, for example the t:N LECs di and do not appear 
in the formalism for NN scattering at the present chiral 
orders. The computational overhead of AD in terms of 
wall time is very small. On a single computational node, 
the calculation of all first- and second-order derivatives of 
the 4450 NN scattering observables with respect to the 
26 LECs at NNLO only takes twice as long as computing 
just the central values. 

It is straightforward, but slightly more costly, to ap¬ 
ply the AD technology to the NCSM diagonalization of 
the nuclear Hamiltonian H for A < 4. If the eigenvalue 
spectrum is non-degenerate, the first-order derivatives of 
the ground-state energy Eq and wave function \ 2 po) with 
respect to the LEC ax are given by [89] 


dEo 


da 


X 



dn 

dax 



d 

dax 


i^o) = E 


(V’il IV’o) 

Eo-Ei 




(26) 

(27) 


Higher-order derivatives are simply obtained by repeated 
differentiation. 

Eor bound-state observables, the computational over¬ 
head in terms of wall time is slightly larger than for two- 
body scattering since we must compute all eigenvalues 
and eigenvectors of R. The calculation of all first and sec¬ 
ond derivatives for all 26 LECs at NNLO for the A = 3 
observables is approximately 20 times slower than just 
calculating the central values. 


G. Uncertainty Quantification 

We employ well-known methods from statistical re¬ 
gression analysis to study the sensitivities and quan¬ 
tify the uncertainties at the optimum x^(a^), see 
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e.g. Dobaczewski et al [79]. The Ncx x covariance 
matrix Cov(a^) defines the permissible variations Aol in 
the LECs that maintain an objective function value such 
that 


x"(a* + Aa)-x'(a*)<T, (28) 


where T is some chosen tolerance. We can assume rather 
small variations Aa, and therefore truncate a Taylor ex¬ 
pansion of the objective function at the second order 


X^(a* + Aa) - x^(a*) -(Aa)^H(Aa), 


where Hij = 




d aid an 


(29) 


observable and compute the linear correlation coef¬ 
ficient between any two observables Oa and Ob- To this 
aim, it is most convenient to use the rotated and inde¬ 
pendent LEG representation x defined above. Each LEG 
Xi is normally distributed with zero mean. Next we use 
a quadratic approximation of the observable Oa-) 

Oa{ol^ + Aa) - Oa{ol^) 

«(Aa)^J^ + i(Aa)^H^(Aa) 

= + ix'^U'^H^Ux 

= x^J^ + ix'^H^x, 


are matrix elements of the Hessian H. This should 
be positive definite. It can be decomposed into H = 
UDU^, where the columns of U are the eigenvectors of 
H and D is a diagonal matrix with the eigenvalues of H. 
Defining x = U^(Aa), Eq. (28) becomes 

1 1 

-x'^Dx = 2 E 

The Noc parameters x can be viewed as “rotated” LEGs. 
They are very convenient since they are independent of 
each other, which simplifies the previous equation and 
gives 

< Ti Vi, (31) 

where Ti is the limit to use when considering only vari¬ 
ations in one parameter and keeping the others fixed. If 
X^(a) follows a chi-squared distribution, then x\D^^^iij2 
will also follow a chi-squared distribution with one de¬ 
gree of freedom, meaning that the la confidence level 
is given by Ti = 1, and Xi ^ A7(0, 2 /D^ 2 ^^^). In prac¬ 
tice, x^(a^) will only be an approximate chi-squared 
distribution, which modifies Ti slightly. Here we set 
Ti = x^(a*)/Adof which corresponds to a rescaling of 
the x^(a^)-function [79], 

Xscaied (a) = X^ (a) • (32) 

The covariance matrix is then given by 

Cov(a*) = = USU^, (33) 

Ndof 

where 5] is the diagonal matrix with the vector of vari¬ 
ances, cr^, of the rotated LEGs, on the diagonal. Since 
Ti only affects Gov with a constant factor, correlations 
remain invariant under changes in Ti. 


where Ja is the Jacobian vector of partial derivatives, 
JA,i = ^^7 Ha is the corresponding Hessian matrix, 
and the tilde notation in the last line indicates the sim¬ 
ilarly rotated Jacobian and Hessian. The corresponding 
statistical expectation value E(-) is given by 

1 ^ 

E[C>^(a)] « C>A(a*) + 

ij (35) 

= OA{a^,) + dmg{HA) 

Einally, we define the covariance of Oa and by 

Cov(A, B) = E[ ( 0 ^(a) - £[0^(0)]) 

X (OB(a)-E[(!)B(a)])] 

~ 1 - 1 - 
~ E® + 2^A,ijXiXj - -HA,iicri) 

ijkl f / 

X {Js.kXk + -T^HB^klXkXl - -HB,kk(^k)\ 

= JJSJb + o Hb)o-2, 

where o denotes the Hadamard product. The statisti¬ 
cal unce rtainty of an observable Oa is then given by 
a A = \/Cov^H 7A). This approximation of the covari¬ 
ance is valid as long as the quadratic approximations (29) 
and (34) are valid and the normalized objective function 
can be assumed to follow a chi-squared distribution. 

Using a linear approximation, the probability distri¬ 
bution for an observable Oa will follow the well-known 
Gaussian form. However, for the quadratic approxima¬ 
tion there is no such analytic expression. Instead, it is 
easy to reconstruct the probability distribution numeri¬ 
cally by using Eq. (34) with a large sample of parameter 
sets. 


1. Error propagation 


III. RESULTS 


Starting from the covariance matrix Gov(a^) we can 
propagate the statistical uncertainties in the LEGs to any 


In this section we discuss our results from the optimiza¬ 
tion of xEET at LO, NLO, and NNLO (Sec. HI A), the 
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subsequent error propagation (Sec. IIIB), as well as an 
expanded discussion on the implications and advantages 
of a simultaneous optimization protocol (Sec. IIIC.) In 
particular we discuss the important consequences of cor¬ 
relations between the LECs in the case of simultaneous 
versus separate optimization strategies. 

A. Optimization 

With all the necessary tools in place we can perform 
the fits to experimental data. For all cases we implicitly 
assume that the LECs are of natural size [39] by choosing 
starting points in this region of the parameter space. We 
did not in any other way force the LECs to be natural. 
A possible problem in multi-parameter optimization is 
the existence of several local minima. At LO, with just 
two parameters, there is only one minimum. However, 

Table III. Comparison of different minima at various chi¬ 
ral orders. AW-LECs are optimized using only NN scat¬ 
tering data (at NNLO, the ttN LECs are fixed). The min¬ 
ima are equally good for A = 2 observables, but differ sig¬ 
nificantly in A = 3 bound-state properties, calculated here 
without a three-body force. The last row corresponds to 
parameters and results (with NN forces only) of the simul¬ 
taneously optimized NNLOsim interaction. The C LECs 
are in units of lO^GeV”^. The scattering y^/iVdof shown 
are for data up to 125 MeV without model errors included. 
^(exp)(3jj^ _ -8.48 MeV. Energies are in MeV. 



Mrip) 


xViVdof 

SfH) 


LOsep 

-0.11 

-0.072 

350 

-2.21 

-11.4 

NLO-1 

+0.81 

+0.69 

14 

-2.17 

-3.03 

NLO-2 

+0.81 

-0.17 

14 

-2.16 

-3.30 

NLO-3 

-0.15 

+0.68 

14 

-2.17 

-2.92 

NLO-4 

-0.15 

-0.17 

14 

-2.16 

-8.22 

NNLO-1 

+0.49 

+0.53 

2.4 

-2.19 

-3.64 

NNLO-2 

+0.49 

-0.17 

2.4 

-2.21 

-3.71 

NNLO-3 

-0.15 

+0.53 

2.4 

-2.19 

-3.23 

NNLO-4 

-0.15 

-0.17 

2.4 

-2.22 

-8.21 

NNLOsim 

-0.15 

-0.17 

1.7 

-2.22 

-8.54 


at NLO we find four local minima. They correspond to 
combinations of two optima in the ^Sq channel and two 
optima in the coupled ^Si —^Di channel. As shown in Ta¬ 
ble III, all four combinations describe scattering data and 
the deuteron properties equally well, thus making them 
indistinguishable from this point of view. Furthermore, 
a similar set of minima exists at NNLO when fitting the 
irN and NN data separately. 

A theoretical argument can provide partial guidance 
in the choice between these parameter sets. The nuclear 
interaction will have an approximate Wigner SU(4) sym¬ 
metry [90] due to the large scattering lengths in the S- 
waves, which implies CiSo ~ This approximate 

constraint rules out the second and third of the four can¬ 
didate NLO and NNLO minima in Table III. Further¬ 
more, we might argue that the fourth minimum (NLO-4 


and NNLO-4, respectively) is the physical one since its 
C LECs most resemble the values obtained at LO. This 
is not a strong justification since LECs are allowed to 
vary between orders. In the end, it does turn out that 
both NLO-4 and NNLO-4 are indeed close to the single 
minimum that exists in the simultaneous NNLO opti¬ 
mization. 

A much more interesting difference between the four 
minima occurs in the few-nucleon sector. It turns 
out that minima 1-3 give significant underbinding of 
the triton. Since the measured ground-state energy 
is -8.48 MeV, these results imply that three-nucleon 
forces, which appear at NNLO, would have to contribute 
5-6 MeV of the missing binding energy. This difference 
is smaller for the NLO-4 and NNLO-4 minima and they 
most likely represent the physical minima. This is also 
more in line with the power-counting arguments that 
the three-nucleon force should be weaker than the two- 
nucleon force, see e.g. Ref [91]. Furthermore, with the 
subsequent addition of the NNN terms at NNLO (as it 
is done in the sequential optimization strategy) it turns 
out that only the NNLO-4 minimum allows to reproduce 
all A = 3 observables within one standard deviation. For 
these reasons, NLO-4 and NNLO-4 define the A/7V-only 
parts of the NLOsep and NNLOsep potentials, respec¬ 
tively. 

The values for the LECs of our optimized potentials at 
LO, NLO, and NNLO are tabulated in the Supplemen¬ 
tal Material [92] together with their estimated statistical 
uncer tainties. T he statistical uncertainty of the Rh LEC, 
i.e. is a measure of how much this particular 

parameter can change while maintaining a good descrip¬ 
tion of the fitted data, as detailed in section IIG. That 
is, the uncertainty for a given LEC represents its max¬ 
imal variation while assuming that all other LECs are 
fixed at the minimum. Note, however, that the LECs 
really cannot be varied independently of each other due 
to mutual correlations. A full error analysis requires a 
complete covariance matrix as we demonstrate below. 

The appearance of NNN diagrams and sub-leading 
terms from the ttN sector does not occur until NNLO 
in our chiral expansion. This implies small differences 
between the separately and simultaneously optimized in¬ 
teractions at lower orders. The deuteron properties are 
included in the optimization of LOsim and NLOsim, but 
not in LOsep and NLOsep. We find that the statistical 

values (not including the model errors) with respect to 
NN scattering data are almost identical for LOsim and 
LOsep, and so are the values of the LECs. The small 
value of CTexp+method foi* the dcutcron binding energy con¬ 
strains the statistical error for Cas^ in LOsim correspond¬ 
ingly. For the contact potential at NLO there are three 
LECs that operate in the deuteron channel, more than 
in any other NN partial wave. The presence of mutual 
correlations cannot be neglected. This explains why the 
individual statistical errors for the LECs in NLOsim and 
NLOsep in this channel are similar and larger than at 
LO. The covariances will also impact the value for the 
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forward error in the deuteron binding energy, discussed 
further in Sec. IIIB. 

We find that the description of the pp scattering data 
is not influenced much by the inclusion of the deuteron 
in the optimization, while the agreement with np data is 
notably worse above 35 MeV. At this order it is mainly 
the Cssi LECs that have changed, see Ref. [92], 

which only affect np scattering. 

As previously mentioned, the yEFT interaction be¬ 
comes significantly more involved at NNLO as NNN and 
sub-leading ttN terms enter at that order. The simulta¬ 
neous optimization of all data listed in Table I leads to 
the construction of the NNLOsim interaction. The conse¬ 
quences of the simultaneous approach are dramatic. First 
of all, we find a single optimum as this strategy elimi¬ 
nates all but one of the local minima that were obtained 
in the sequential optimization. Moreover, a possible con¬ 
cern turns out to be unwarranted: an improved overall 
description of scattering data does not detoriate the de¬ 
scription of different subsets. In fact, the result is quite 
the opposite. With the simultaneous-optimization strat¬ 
egy we find that the description of the pp scattering data 
is actually significantly improved. For scattering ener¬ 
gies Tiab < 290 MeV the statistical xfpp)/^doi = 9.1 for 
NNLOsim compared to xfpp)/^dof = 26 for NNLOsep, 
not including the model error. At the same time, the 
for np scattering and irN scattering are similar for the 
two potentials. Measured np scattering cross sections 
are characterized by larger uncertainties and it is there¬ 
fore not surprising that this data remains well described. 
However, it is noteworthy that the NNLOsim potential 
reaches a better description of the NN data while main¬ 
taining a description of the ttN data that is comparable 
to the one of NNLOsep. Keep in mind that NNLOsep is 
separately optimized to the irN scattering data. In the 
simultaneous optimization protocol we are effectively in¬ 
troducing additional constraints on the q LECs via the 
NN data set. One could be concerned that the short- 
range NN physics would impact and worsen the descrip¬ 
tion of the long-range pion physics. It is not unlikely 
that we would have seen such unphysical effects if the 
ttN database would have been more comprehensive. The 
existing irN data does not constrain all directions in the 
ttN leg parameter space, which allows for large varia¬ 
tions in the parameter values and a better description of 
the pp scattering data with NNLOsim. The y^/A^dof for 
NN and irN scattering up to different are presented 
in Fig. 5. 

The predominant advantage of the simultaneous op¬ 
timization is the correct treatment of correlations. Al¬ 
though the uncertainties of the LECs presented in the 
Supplemental Material [92] are similar for NNLOsep 
and NNLOsim, the propagated statistical errors of ob¬ 
servables can be several orders of magnitude larger for 
NNLOsep due to missing correlations, see Sec. IIIB. To 
visualize the correlations between all LECs, we plot the 
linear correlation matrix in Fig. 6 . The linear correla¬ 
tion between two LECs, or any observables A and H, 





Figure 5. (a) Cumulative y^/Ndof for NN scattering data 

including the model error (see Sec. HE). Note that the am¬ 
plitude of the model error is chosen so that y^/Ndof = 1 when 
all data up to Tiab = 290 MeV is included. (b,c) Cumulative 
y^/Vdof without the model error for NN and ttN scattering 
data, respectively. 


indicates their linear relationship and is defined as the 
normalized covariance, Coy{A, This quan¬ 

tity assumes values between —1 (fully anti-correlated) 
and +1 (fully correlated). A positive (negative) value for 
the correlation indicates that a larger value for A most 
likely requires a larger (smaller) value for B. The corre¬ 
lation coefficients between LECs that belong to different 
objective functions are zero. For NNLOsep this implies 
that the correlation matrix is block-diagonal in terms of 
the ttV, a/TV, and NNN sectors. For NNLOsim, however, 
such inter-block correlations are revealed. In addition, we 
observe an increase of the correlations within each group. 
This can be traced to the fact that the ttTV LECs, ci, C 3 
and C 4 , occur in the description of TVTV-, ttTV-, and NNN- 
data. The failure to capture these correlations within 
the sequential optimization approach, such as with the 
NNLOsep potential, will induce very large propagated 
statistical errors. In conclusion, simultaneous optimiza¬ 
tion is key for a realistic forward propagation of para¬ 
metric uncertainties. 












14 


OOOOtHt-h OtH^hcm 

CO COCOCO cocorJc^i^i^CL, 

,H tHtHtH COCOl^COrHPOCO 

<0 'O'O O'OOOOOOO 


CJ CJ 


io i;o 

CM CO ^ tH iH 

O O CJ 0) 


Cl50nn 

0000000000 









ClSOnp 


■ 

✓ 


0000000 









ClSOpp 


■ 

■ 

\ 

0000000 









Ciso 


■ 

■ 

■ 

OIOIOIOOPO 









Cssi 








0000 









Cssi 





■ 

■ 

00000 









Cei 





■ 


■ 

0000 









Cspo 








■ 

000 









CiPi 









■ 

00 









Cspi 










■ 

0 









C 3 P 2 








pi 

iM| 


m 

I 








CD 

A 

fA 










■ 

3 







CE 













m 

1 






Cl 









} 

VI 

iH 


n 

■ 


000 


C2 














■ 

■ 

000 

\ 

C3 
















■ 

000 

C4 

















■ 

00 

ei5 


















■ 

0 

ei6 



1 









TT 

]vj 

0 




■ 


e a, 
g e 
o o 
CO CO 


CO CO CO CO 


a, a, a, a, 


'0*0*0 0*0 o o o o o o 


Q ^ 

O O CJ 



/ 


/ 


0 

\^\\ 

00 

// 

\oo\ 

■ 

/ 

z 

✓ 

z 

0 

\^\\ 

00 

// 

\oo\ 

■ 

■ 


✓ 

z 

0 

\'^\\ 

00 


\oo\ 

■ 

■ 

■ 

✓ 


ooz|zMo|o\ 

■ 

■ 

■ 

■ 

✓ 

0 

0 


00 

00 

✓ /’ 

\ 00 '^ 

\oo\ 


|OOOOOc:iOOO^OO 




■ 

■ 

■ 

■ 


■ 

✓ 

z 

z 


0 

\ 

\ 

ZO^ 

z 





□ 



■ 

■ 

/ 

z 


0 



/oo 

/ 











z 

^0 

\Nzroio 

z 













0 

\ 

\ 

zoo 

z 













|o 

0 

0000 


1000^00 

l/\o^\ 
\\0(p\ 

moo 


I 


1.0 

0.8 

0.6 

0.4 

0.2 

0.0 


03 


<D 

0.2 I 


-0.4 

- 0.6 

- 0.8 

- 1.0 


Figure 6. Graphical representation of the linear correlation matrix for NNLOsep (left) and NNLOsim (right) including selected 
LECs. The separately optimized NNLOsep potential does not probe the statistical correlation between LECs entering different 
optimization stages. It is striking that there are almost no correlations for the NNLOsep potential, while for the NNLOsim- 
potential the situation is quite the opposite. 


B. Error propagation 

Statistical errors and covariances between computed 
observables are calculated under the assumption that 
each observable depends quadratically on the LECs in 
the vicinity of the minimum, see Eq. (34). Our estimate 
of the statistical uncertainty, a a, of an observable, Oa, 
rests on this assumption, which also explains why we 
have asymmetric error bars. We have performed exten¬ 
sive Monte Carlo samplings to verify the validity and ne¬ 
cessity of using the second-order approximation. A linear 
truncation is more common. In particular, we compare 
the probability density function for various observables 
obtained from: (i) Monte Carlo samplings of the multi¬ 
variate Gaussian spanned by the covariance matrix, (ii) 
the quadratic approximation, and (iii) the linear approx¬ 
imation of Eq. (34). The Monte Carlo calculations use 
10 ^ sets of normally distributed LEC vectors. 

The probability distributions for the scattering lengths 
app and for the potentials NNLOsep and NNLOsim 
are shown in Fig. 7. Note that these results are predic¬ 
tions since the scattering lengths are not included in the 
objective function at NNLO. The statistical errors for 
and obtained in the Monte Carlo calculations with 
the NNLOsim potential are small and well reproduced 
already by the corresponding linear approximation, as 
expected. With NNLOsep, the errors are much larger 
and require at least a quadratic approximation for the 
forward error. The uncertainties of the ERE parameters 


differ quite a lot between these two potentials. It is im¬ 
portant to remember that for the NNLOsim potential, 
all LECs are constrained by irN, NN as well as NNN 
data. Hence, in the error analysis, the LECs that fulfill 
^scaled® ^ ^dof will provide a reasonable description of 
most scattering data. The ttN LECs for NNLOsep on the 
other hand, are constrained only by the TrWdata and the 
missing statistical correlations allow for wide permissible 
ranges for the NN scattering lengths. 

It is possible to explore correlations between any pair 
of observables by looking at joint probability distribu¬ 
tions. As an example, we plot the statistical distribution 
of binding energies of ^He and corresponding radii of the 
deuteron for the NNLO potentials in Fig. 8 . The con¬ 
tour lines indicate the regions that encompass 68 % (la) 
and 95% (2a) of the probability density. It is remarkable 
that the quadratic approximation (dashed lines) repro¬ 
duces even the fine details of the full calculation (solid 
lines) for the NNLOsim interaction. Again, the magni¬ 
tude of variations is strikingly large for NNLOsep, but 
the quadratic approximations does rather well in repro¬ 
ducing them. In particular, we see a large improvement 
when going from a linear (dotted lines) to a quadratic de¬ 
pendence on the LECs. This even captures the departure 
from the standard first-order ellipse. 

We present final results for bound-state observables in 
few-body systems (A = 2 — 4) as well as ERE parame¬ 
ters in Table IV for the LO, NLO and NNLO potentials. 
Observables that were part of the respective objective 
function are indicated by a white background, while en- 
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Figure 7. Histograms (filled green area) for the sampled 
probability distribution of the nn (b,d) and pp (a,c) scatter¬ 
ing lengths (including Coulomb) using the NNLO potentials: 
NNLOsep (a,b) and NNLOsim (c,d). The dashed (solid) lines 
show error estimates from the sample assuming that the scat¬ 
tering length depends linearly (quadratically) on the htting 
parameters. The hnal theory result (red square) from Eq. (34) 
agrees well with the sampled distribution. 


tries with grey background are predictions. Note that 
the errors that are given in this table do not include a 
model error from the xEFT truncation, only the propa¬ 
gated statistical uncertainties as described in Sec. IIG. It 
is therefore difficult to make strong conclusions regarding 
the order-by-order convergence but we certainly observe 
improved predictions when going to higher orders. Note 
that LO results in general are characterized by small sta¬ 
tistical uncertainties since very little freedom is allowed 
with just two parameters. At NNLO we observe large 
statistical errors in the predictions following the sequen¬ 
tial approach, e.g., with NNLOsep the statistical error 
for E(^He) is more than 10 MeV. 

Energies and radii of few-nucleon systems are well 
reproduced by NNLOsim as shown in Table IV, with 
the deuteron radius being the possible exception. This 
can be traced back to omitted relativistic effects. Eor 
the deuteron, Ar^ has been estimated to be of the size 
0.013 fm^ [93] and 0.016 fm^ [94]. 

We have also extracted correlations between other ob¬ 
servables in the few-nucleon sector. As expected, for both 



Ei^Ue) (MeV) E{^Ue) (MeV) 


Figure 8. Joint statistical probability distribution for E(^He) 
and rpt_p(^H) for (a) NNLOsim and (b) NNLOsep obtained 
in a Monte Carlo sampling (Agampie = 10^) as described in 
the text. Contour lines for this distribution are shown as 
black, solid lines, while blue dotted (red dashed) contours are 
obtained assuming a linear (quadratic) dependence on the 
LECs for the observables. 


NNLOsep and NNLOsim there exist a significant correla¬ 
tion between the D-state probability and the quadrupole 
moment of the deuteron. More interestingly, at the 
present optima the triton ^^-decay half-life does not cor¬ 
relate strongly with any other bound-state observable in 
Table IV. This corroborates the importance of using this 
observable to constrain nuclear forces, as was done al¬ 
ready in Ref. [77]. 

Total uncertainties (statistical plus estimated model 
error from the xEFT truncation) are shown for scattering 
observables in Fig. 9. The statistical errors are typically 
very small compared to the model error for the LOsim, 
NLOsim and NNLOsim potentials. Fig. 9(b,c,d,e). Clear 
signatures of an order-by-order convergence are seen in 
the NN scattering observables as illustrated by the np 
total cross section and the differential cross section that 
are shown in Fig. 9(b,c). The same convergence is not 
seen when using the sequentially optimized potentials as 
illustrated in Fig. 9(a). In this case, the statistical errors 
are of the same order of magnitude as the model errors, 
and the NNLO error band is even wider than the NLO 
band. Note that irN scattering is only described with 
the NNLOsep and NNLOsim interactions. 


C. Optimization protocol 

We have demonstrated that the statistical uncertain¬ 
ties of xEET, if all correlations are accounted for, will 
induce rather small errors in the predictions of observ¬ 
ables. This reflects the fact that most of the few-nucleon 
data is precise and diverse enough to constrain a statisti¬ 
cally meaningful x^ET description of the nuclear inter¬ 
action. Also note that we only included experimentally 
observable data in the objective function. 

The existence of strong correlations between the LECs 
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Table IV. Statistical uncertainties propagated from the NN, NNN, and ttN LECs to the ground-state energies (in MeV) and 
radii (in fm) for ^ < 4 nuclei, the deuteron D-state probability T)(^H) (in percent) and quadrupole moment Q(^H) (in fm^) and 
effective range observables for the ^Sq channel (in fm). Gray background indicates that the corresponding result is a prediction. 
Asymmetrical errors are due to the quadratic dependence of the observables on the LECs. The error bars on the experimental 
values for bound-state observables include both experimental and method uncertainties as detailed in Table II. 



LOsep 

NLOsep 

NNLOsep 

LOsim 

NLOsim NNLOsim 

Exp. 

Ref. 

eCh) 

-2.211(15) 

-2-163|l^') 

9 9 (+ 12 ) 
^•^(-25) 

-2.223 

-2.224]!g] -2.224|!°] 

-2.225 

Table II 

eCb) 

-11.40(4) 

-8.220 l^^j 

-8 

-11.43 

-8.268 !“] -8.482 !“] 

-8.482(28) 

Table II 

S(®He) 

-10.39(4) 

-7.474|+f4 

7 7(+30) 

‘ (-62) 

-10.43 

-7.528|!“] -7.717|!^I] 

-7.718(19) 

Table II 

S(^He) 

-40.27(13) 

-27.56[!J^3> 

9o(+8) 


-40.38(1) 

-27.44]!“] -28.24]!f^) 

-28.30(11) 

Table II 

rpt-p("H) 

+1.916(5) 


_l_l gyNsr) 
-hi.a ( (_52) 

-hl.912 

+ 1 . 972 ]!°] +1.966]!°] 

+1.976(1) 

Table II 

rpt-p(®H) 

+1.293(2) 

+1.596(3) 

+1.58+^^] 

-hi.292 

+1.614]!3] +1.581(2) 

+1.587(41) 

Table II 

rpt-p("He) 

+1.370(2) 

+1.778[+®j 

+1.76+“) 

-hl.368 

+1.791(3) +1.761(2) 

+1.766(13) 

Table II 

rpt-p("He) 

+1.081(1) 

+1.459(4) 

-hi 44*^^^^^ 


-hl.080 

+1.482(3) +1.445(3) 

+1.455(7) 

Table II 

E\eB) 

- 

- 

+0.685 +“ 

- 

+0.6848(11) 

+0.6848(11) 

Table II 

dCe.) 

+7.794(17) 

+2.942[+|f, 

+3.9h;«) 


-h7.807 

+2.876]!®^] +3.381]!“] 

- 


QfH) 

+0.3035(7) 

+0.2602[+^®) 

+0.2701!“] 

-h0.3030 

+0.2589]!“] +0.2623(8) 

+0.270(11) 

Table II 

'^nn 

-26.04(5) 

-18-95|Tv] 

_19(+7) 

-'■^(-24) 


-26.04(8) 

-18.95]!°°] -19.28]!^^g] 

-18.95(40) 

[18] 

ujnp 

-25.58(5) 

90 qyC+lS) 
Z9.a/(_ig) 

_24^+iO 

^^(-44) 


-25.58(8) 

-23.60]!]°] -23.83(11) 

-23.71 

[95] 

C 

Upp 

-7.579(4) 

-7.799[+^) 

7 qNio) 

< •0(_24) 


-7.579(6) 

- 7 . 799 ]!]] -7.811(1) 

-7.820(3) 

[62] 

' nn 

+1.697 

+2.752(7) 

+2.85[!“] 


+1.697(1) 

+ 2 . 752 ]!]] +2.793(14) 

+2.75(11) 

[18] 

' np 

+1.700 

+2.650[+4j 

+2.74!“] 


+1.700(1) 

+2.648(3) +2.686(2) 

+2.750(62) 

[95] 

pC 

' PP 

+1.812 

+2.704(3) 

+2.8l|!f„] 


+1.812(1) 

+2.704(3) +2.758(2) 

+2.790(14) 

[62] 


require a complete determination of the corresponding 
covariance matrix, not just the diagonal entries. For this, 
it is necessary to employ the so-called simultaneous opti¬ 
mization protocol. To further demonstrate this point, we 
carried out error propagations with NNLOsim while ne¬ 
glecting the off-diagonal correlations between the LECs. 
The statistical uncertainty of the binding energy in ^He 
grew with a factor ^ 90 compared with the fully informed 
model. Neglecting the statistical correlations will also 
obscure the desired convergence pattern of xEFT. In¬ 
deed, for the separately optimized potentials there were 
no signs of convergence in the description of, e.g., np 
scattering data. 

If the experimental database of irN scattering cross 
sections would be complete, then it would be possible to 
separately constrain, with zero variances, the correspond¬ 
ing LECs. Only this scenario would render it unnecessary 
to include the irN scattering data in the simultaneous 
objective function. Implicitly, this scenario also assumes 
a perfect theory, i.e. that the employed xEET can ac¬ 
count for the dynamics of pionic interactions. Of course, 
reality lies somewhere in between, and a simultaneous 
optimization approach is preferable in the present situ¬ 
ation. There exists ongoing efforts where the ttN sector 
of xEET is extrapolated and fitted separately in the un¬ 
physical kinematical region where it exhibits a stronger 
curvature with respect to the data [96] . 

Overall, the importance of applying simultaneous opti¬ 
mization is most prominent at higher chiral orders, since 


the sub-leading ttN LECs enter first at NNLO. In fact, 
the separately optimized NNLOsep potential contains a 
large systematic uncertainty by construction. We find 
that the scaling factor for the NN scattering model error, 
Cnn^ decreases from 1.6 to 1.0 mb^/^ when going from 
NNLOsep to the simultaneously optimized NNLOsim. 
This implies that the separate, or sequential, optimiza¬ 
tion protocol introduces additional artificial systematic 
errors not due to the chiral expansion but due to incor¬ 
rectly fitted LECs. This scenario is avoided in a simulta¬ 
neous optimization. The scaling factor for the irN scat¬ 
tering model error, remains at 3.6 mb^/^ for both 

NNLOsep and NNLOsim. 

The size of the model error is determined such that 
the overall scattering x^/^dof is unity, which means that 
it depends on the observables entering the optimiza¬ 
tion. We can explore the stability of our approach by 
re-optimizing NNLOsim with respect to different trun¬ 
cations of the input AWscattering data. To this end 
we adjust the allowed between 125-290 MeV in six 
steps. It turns out that our procedure for extracting the 
model error is very stable. The resulting normalization 
constants Cnn vary between 1.0 mb^/^ and 1.3 mb^/^ as 
shown in Eig. 10(a). 

To see the corresponding effect on predicted observ¬ 
ables we consider the np total cross section at laboratory 
scattering energy Tiab = 300 MeV. The model errors vary 
between 4.8 mb and 6.1 mb, and the calculated cross sec¬ 
tions vary between 36.5mb and 42.7mb, see Eig. 10(b). 
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Figure 9. Comparison between selected NN and ttN experimental data sets and theoretical calculations for chiral interactions 
at LO, NLO and NNLO. The bands indicate the total errors (statistical plus model errors), (a) np total cross section for the 
sequentially optimized interactions with no clear signature of convergence with increasing chiral order. All other results are for 
the simultaneously optimized interactions: LOsim, NLOsim and NNLOsim. (b) np total cross section; (c) np differential cross 
section; (d) ttN charge-exchange, differential cross section; (e) ttN elastic, differential cross section. 
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Figure 10. Predictions for the different re-optimizations of 
NNLOsim. On the x-axis is the maximum Tiab for the NN 
scattering data used in the optimization, (a) Model error 
amplitude (20) re-optimized so that x^/^dof = 1 for the re¬ 
spective data subset, (b) Model prediction for the np total 
cross section at Tiab = 300 MeV with error bars representing 
statistical and model errors for the different re-optimizations. 


The measured value is 34.563(174) mb [21, 97]. We uote 
that the size of the estimated model error is compara¬ 
ble with the variatiou iu the predictious due to chaugiug 

q^max 
^lab • 

Throughout the aualysis, the model error for scatter- 
iug observables was assumed to scale with momeutum p 
accordiug to Eq. (20). However, the soft scale Q iu yEFT 


is set by max{p, aud it cau be argued that the model 
error should be implemeuted as 


~ (amp) _ ^ 
^ model,X ~ 


/ max{p, 

[ W 




(37) 


It turus out that resolviug these two momeutum scales 
has a small impact ou the estimated model errors. As 
au illustratiou, the predictious of the ^He biudiug euergy 
chauges by just ^ 20keV (less thau 0.1%). Iu fact, this 
effect is much smaller thau the impact of chaugiug the 
cutoff iu the experimeutal NN scatteriug database. 


IV. EXTENDED ANALYSIS OF SYSTEMATIC 
UNCERTAINTIES 

Iu unclear physics the theoretical uncertainties very 
often dominate over the experimental ones. In partic¬ 
ular, this is true for the systematic error. Therefore, 
it is crucial to establish a credible program for assess¬ 
ing the error budget of any prediction or analysis of 
experimental information. Thus, we focus our atten¬ 
tion on the convergence and missing physics in yEFT. 
In particular, we discuss consequences for predictions of 
bound-state observables in heavier nuclei such as ^He 
and It would be valuable to estimate the systematic 
uncertainty of predicted bound-state observables — due 
to the momentum-dependent yEFT uncertainty (Tmodei 
in Eq. (20). However, the explicit momentum depen¬ 
dence is integrated over when solving the non-relativistic 
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Schodinger equation. Thus, a clear connection to the 
momentum-expansion is lost. 

As demonstrated already in Fig. 10, the variations in 
model predictions obtained from different truncations of 
the input data (including only NN scattering data with 
^lab ^ ^ good first approximation of the ex¬ 

pected model uncertainty. To get a more complete pic¬ 
ture of the systematic uncertainty, we now also vary the 
regulator cut-off parameter A in the range 450 — 600 MeV 
in steps of 25 MeV. For each combination of and 

A we perform a simultaneous optimization of the LECs, 
which results in a family of 42 NNLO interactions - i.e., 
42 sets of LECs that each comes with statistical uncer¬ 
tainties. It is clear from Table V that the statistical 
uncertainties of the LECs are smaller than the overall 
shifts induced by varying and the cutoff A. All sets 
of LECs at LO, NLO, NNLO that were obtained in this 
work are listed in the Supplemental Material [92]. Eur- 
thermore, each set is accompanied by its own covariance 
matrix, also available for download. In the following dis¬ 
cussion we use this family of potentials to estimate the 
systematic uncertainty. 

Eirst, we would like to emphasize that all sets of si¬ 
multaneously optimized LECs provide an almost equally 
good description of all A < 4 data. Some of the ttN 
LECs display large variations, but the x^/A'dof (without 
model error) for the ttN data is within 2.28(4) for all of 
these potentials. The sub-leading ttN LECs become more 
positive when NN scattering data at higher energies is 
included, and ci in particular carries a larger (relative) 
statistical uncertainty than the others. It is noteworthy 
that for a given and up to la precision, the irN 

LECs exhibit A-independence. The NNN LECs, cd and 
ce, tend to depend less on at larger values of A. 

However, they always remain natural. It is also interest¬ 
ing to note that the tensor contact, Cei, is insensitive to 
A-variations but strongly dependent on the cut. It 

was shown in Eig. 6 that Ce^ and C 4 correlate strongly. 
This effect can already be expected from the structure of 
the underlying expression for the NNLO interaction. 

To gauge the magnitude of model variations in heavier 
nuclei we computed the binding energies of ^He and 
using the previously mentioned family of 42 NNLO po¬ 
tentials. The resulting binding energies for ^He and ^^O, 
computed in the NCSM and CC, respectively, are shown 
in Eig. II. The NCSM calculations were carried out in 
a HO model space with AViax = 20 and Tiuj = 36 MeV. 
The CC calculations were carried out in the so-called 
A—CCSD(T) approximation [7] in 15 major oscillator 
shells with hw = 22 MeV. The largest energy difference 
when going from 13 to 15 oscillator shells was 3.6 MeV 
(observed for A = 600MeV). Eor our purposes, this pro¬ 
vides well-enough converged results. The NNN force was 
truncated at the normal-ordered two-body level in the 
Hartree-Eock basis. 

The E(^He) predictions vary within a ^ 2 MeV range. 
Eor F^(^^O) this variation increases dramatically to ^ 
35 MeV. Irrespective of the discrepancy with the mea- 


Table V. Ranges of LEG values and maximum statistical un¬ 
certainties among all 42 simultaneously optimized NNLO po¬ 
tentials constructed in this work, see Sec. IV The first two 
columns show the global variation of the LEG values, in terms 
of minimum and maximum values, due to changes in A and 
The third column shows the maximum statistical un¬ 
certainty of each LEG, which almost exclusively come from 
the A = 450 MeV and = 125 MeV NNLOsim potential. 
For a given LEG, the statistical uncertainty is rather similar 
for different potentials. Ci are in units of lO^GeV”^, Ci in 
units of 10^ GeV“^, cd and ce are dimensionless while Ci, di 
and ei are in units of GeV“^, GeV“^ and GeV“^, respectively. 


LEG 

range 


max(cr) 

/Xf(np) 

-0.I5I9 


-0.1464 

±0.0020 

Mpp) 

-0.I5I2 


-0.1454 

±0.0020 

^{nn) 

-0.I5I8 


-0.1463 

±0.0021 


2.4188 


2.5476 

±0.0511 


-0.1807 


-0.1348 

±0.0032 


0.5037 


0.7396 

±0.0521 

Cei 

0.2792 


0.6574 

±0.0253 

Cspq 

0.9924 


1.6343 

±0.0428 

Cip, 

0.0618 


0.6635 

±0.0438 

Csp^ 

-0.9666 


-0.4724 

±0.0416 

Cap, 

-0.7941 


-0.6324 

±0.0327 

CD 

-0.5944 


0.8348 

±0.0833 

Ce 

-2.4019 


-0.0893 

±0.2282 

Cl 

-0.8329 


0.2784 

±0.3043 

C2 

2.7946 


5.3258 

±1.0754 

C3 

-4.3601 


-3.4474 

±0.1506 

C4 

1.8999 


4.2353 

±0.2179 

di +0^2 

4.4636 


5.4505 

±0.1378 

ds 

-4.8549 


-4.4583 

±0.2302 

d5 

-0.2992 


0.0233 

±0.1407 

di4 —di5 

-10.3220 


-9.6902 

±0.2820 

ei4 

-0.3700 


0.9569 

±0.9079 

ei5 

-11.9223 


-9.1307 

±2.4962 

C16 

-0.6847 


7.4463 

±4.2436 

cir 

0.9322 


1.4986 

±1.8143 

CIS 

-2.5068 


8.3777 

±1.9022 


sured value, the spread of the central values indicates 
the presence of a surprisingly large systematic error when 
extrapolating to heavier systems. 

The statistical uncertainties remain small: tens of keV 
for ^He and a few hundred keV for These uncer¬ 
tainties are obtained from the quadratic approximation 
with the computed Jacobian and Hessian for ^He, while 
a brute-force Monte Carlo simulation with 2.5 x 10 ^ CC 
calculations was performed for This massive set of 
CC calculations employed the doubles approximations in 
9 major oscillator shells. We conclude that the statistical 
uncertainties of the predictions for E(^He) and E(^^O) 
at NNLO are much smaller than the variations due to 
changing A or However, this is only true for simul¬ 

taneously optimized potentials. Eor the separately opti¬ 
mized NNLO potential (NNLOsep) the statistical uncer¬ 
tainty of the E(^He) prediction is five times larger than 
the observed variations due to changing A and 
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Figure 11. Binding-energy predictions for (a) ^He and (b) 
with the different re-optimizations of NNLOsim. On the x- 
axis is the employed cutoff A. Vertically aligned red markers 
correspond to different for the NN scattering data used 

in the optimization. The experimental binding energies are 
F^(^He) —28.30 MeV, represented by a grey band in panel 

(a), and E(^^0) ^ —127.6MeV [98]. Statistical error bars 
on the theoretical results are smaller than the marker size on 
this energy scale. 
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V. OUTLOOK 

The extended analysis of systematic uncertainties pre¬ 
sented above suggests that large fluctuations are induced 
in heavier nuclei (see Fig. 11). Furthermore, while pre¬ 
dictions for ^He are accurate over a rather wide range of 
regulator parameters, the binding energy for turns 
out to be underestimated for the entire range used in this 
study. In fact, there is no overlap between the theoreti¬ 
cal predictions and the experimental results, even though 
the former ones have large error bars. 

Based on our findings we recommend that continued 
efforts towards an ab initio framework based on xEFT 
should involve additional work in, at least, three different 
directions: 

1. Explore the alternative strategy of informing the 
model about low-energy many-body observables. 

2. Diversify and extend the statistical analysis and 
perform a sensitivity analysis of input data. 

3. Continue efforts towards higher orders of the chiral 
expansion, and possibly revisit the power counting. 

Let us comment briefly on these research directions. The 
poor many-body scaling observed in Eig. 11 was prag¬ 
matically accounted for in the construction of the so- 
called NNLOsat potential presented in Ref. [35] where 
also heavier nuclei were included in the fit. The accu¬ 
racy of many-body predictions was shown to be much 
improved, but the uncertainty analysis is much more dif¬ 
ficult within such a strategy. 

Secondly, to get a handle on possible bias in the sta¬ 
tistical analysis due to the choice of statistical tech¬ 


nique, it is important to apply different types of opti¬ 
mization and uncertainty quantification methods. Vari¬ 
ous choices exist, such as e.g. Lagrange multiplier anal¬ 
ysis [99], Bayesian methods [100], or Gaussian process 
modeling [101, 102]. In general, stochastic modeling with 
Monte Carlo simulations offer a straightforward and ver¬ 
satile approach. This tool is also indispensable for com¬ 
puting the posterior probabilities in Bayesian inference. 
The Monte Carlo results for A < 4 observables that were 
presented in this work consist of 10^ sampling points 
over a multivariate Gaussian parameter space. With our 
current implementation, the computational cost for sam¬ 
pling all A < 4 observables presented in this work is very 
low — less than 8000 CPU hours. As such, the present 
work shows great promise also for future stochastic ap¬ 
plications. 

Eurthermore, the computational framework that we 
have presented here, and our present implementation, is 
not limited to any particular type of regulator function 
or flavour of chiral expansion. Moreover, the handling of 
a larger number of LECs, as would be the consequence 
of working at a higher chiral order, should be relatively 
straightforward and we don’t foresee any computational 
bottlenecks. 

Einally, the magnitude of the systematic uncertainties 
that were observed in this work suggest the need to fur¬ 
ther explore and improve the theoretical underpinnings 
of the chiral expansion of the nuclear interaction. 
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Appendix: Supplemental material for 
“Uncertainty analysis and order-by-order 
optimization of chiral nuclear interactions” 

Here we present all optimized central values and sta¬ 
tistical uncertainties for the LECs in the ttN, NN and 
the NNN sector for all potentials referred to in the main 
text. The q, di and are in units of GeV“^, GeV“^ and 
GeV“^ respectively. Ci are in units of 10^ GeV“^, Ci in 
units of 10^ GeV“^ and cd and ce are dimensionless. 

The covariance matrices for all potentials and numeri¬ 
cal tables with LECs can be downloaded as a tarball from 
the supplemental material page at the published version 
of the manuscript. They are also available upon request 
from the authors. The ordering of the LECs follow the 
ordering in the tables below. 
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Table VL LOsep and LOsim (A = 500,Tmax = 290) 


LEG 

LOsep 

LOsim 


-0.1076841(50) 

-0.1076845(80) 


-0.07172(11) 

-0.0718086(27) 


Table VIL NLOsep and NLOsim (A = 

500,Tn,a>c = 290) 

LEG 

NLOsep 

NLOsim 


-0.150533(96) 

-0.150623(79) 

Mpp) 

-0.14893(11) 

-0.14891(11) 

^(nn) 

-0.14992(27) 

-0.14991(27) 


+1.6926(82) 

+1.6935(83) 


-0.1742(20) 

-0.1843(16) 

Casi 

-0.408(23) 

-0.218(14) 

Cex 

+0.238(14) 

+0.263(16) 


+1.3085(86) 

+1.2998(85) 

Cl Pi 

+0.849(47) 

+1.025(59) 

C3P1 

-0.3409(98) 

-0.336(10) 

Csp^ 

-0.2011(15) 

-0.2029(15) 


Table VIIL NNLOsep and NNLOsim (A = 500,Tmax = 290) 


LEG 

NNLOsep 

NNLOsim 


-0.15387(10) 

-0.1474(20) 

Mpp) 

-0.152935(72) 

-0.1465(20) 

^(nn) 

-0.15354(43) 

-0.1471(20) 


+2.7442(19) 

+2.548(47) 


-0.1671(10) 

-0.1687(21) 


+0.8738(64) 

+0.705(47) 

Cei 

+0.6899(67) 

+0.597(11) 

C^Po 

+1.2782(66) 

+1.161(31) 

Cip^ 

+0.521(12) 

+0.520(33) 

Csp^ 

-0.9378(69) 

-0.955(31) 

Csp^ 

-0.68645(76) 

-0.658(30) 

CD 

-0.581(28) 

-0.325(51) 

Ce 

-0.6666(99) 

-0.521(17) 

Cl 

-0.69(50) 

+0.22(30) 

C2 

+3.0(14) 

+5.1(10) 

C3 

-4.12(32) 

-3.56(13) 

C4 

+5.35(81) 

+3.933(85) 

di+d2 

+6.22(44) 

+5.320(94) 

ds 

-5.31(30) 

-4.83(22) 

ds 

-0.46(18) 

-0.24(14) 

di4 — di5 

-11.00(42) 

-10.23(27) 

ei4 

-0.63(95) 

-0.26(89) 

ei5 

-7.7(26) 

-9.3(24) 

ei6 

+5.9(49) 

-0.0(41) 

cir 

+2.1(18) 

+1.5(18) 

Cl8 

-8.1(42) 

-1.2(16) 
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Table IX. NNLOsim (A = 450) 


LEG 

TiX = 125 

TiX = 158 

Ta"" = 191 

Ta"" = 224 

Ta"" = 257 

Ta"" = 290 


-0.1519(20) 

-0.1512(20) 

-0.1510(20) 

-0.1501(20) 

-0.1499(20) 

-0.1496(20) 

Mpp) 

-0.1512(20) 

-0.1504(20) 

-0.1502(20) 

-0.1493(20) 

-0.1491(20) 

-0.1488(20) 

P^(nn) 

-0.1518(21) 

-0.1510(21) 

-0.1508(21) 

-0.1498(20) 

-0.1496(20) 

-0.1493(20) 


+2.498(51) 

+2.480(50) 

+2.477(49) 

+2.511(48) 

+2.518(48) 

+2.527(48) 

C^s, 

-0.1743(23) 

-0.1780(22) 

-0.1785(22) 

-0.1801(21) 

-0.1806(21) 

-0.1807(21) 

C^s, 

+0.676(52) 

+0.682(51) 

+0.677(50) 

+0.722(49) 

+0.730(48) 

+0.740(48) 

Cei 

+0.313(25) 

+0.475(20) 

+0.478(18) 

+0.600(15) 

+0.629(14) 

+0.657(13) 

C^Po 

+0.992(37) 

+1.122(34) 

+1.111(33) 

+1.164(32) 

+1.171(32) 

+1.183(31) 

Cl Pi 

+0.064(44) 

+0.333(40) 

+0.355(39) 

+0.551(36) 

+0.610(36) 

+0.664(35) 

Cspi 

-0.893(34) 

-0.940(33) 

-0.943(33) 

-0.960(32) 

-0.961(32) 

-0.967(31) 

Csp^ 

-0.794(33) 

-0.706(32) 

-0.708(32) 

-0.663(31) 

-0.651(31) 

-0.639(31) 

CD 

+0.359(83) 

-0.109(68) 

-0.089(63) 

-0.432(56) 

-0.508(54) 

-0.594(52) 

Ce 

-0.089(29) 

-0.281(28) 

-0.276(26) 

-0.443(24) 

-0.483(23) 

-0.528(22) 

Cl 

-0.83(30) 

-0.53(30) 

-0.50(30) 

-0.21(30) 

-0.13(30) 

-0.05(30) 

C2 

+2.8(11) 

+3.1(11) 

+3.2(11) 

+3.8(11) 

+4.0(11) 

+4.2(11) 

C3 

-4.36(15) 

-3.82(14) 

-3.82(14) 

-3.57(14) 

-3.51(14) 

-3.45(14) 

C4 

+1.90(22) 

+3.02(16) 

+2.95(15) 

+3.84(12) 

+4.02(11) 

+4.235(98) 

di+d2 

+4.47(14) 

+4.91(12) 

+4.88(11) 

+5.28(10) 

+5.36(10) 

+5.450(97) 

ds 

-4.49(23) 

-4.63(23) 

-4.61(23) 

-4.78(22) 

-4.82(22) 

-4.85(22) 

d5 

+0.02(14) 

-0.14(14) 

-0.13(14) 

-0.25(14) 

-0.27(14) 

-0.30(14) 

di4 — di5 

-9.71(28) 

-9.98(28) 

-9.95(27) 

-10.21(27) 

-10.26(27) 

-10.32(27) 

ei4 

+0.96(91) 

+0.33(90) 

+0.35(90) 

-0.14(89) 

-0.25(89) 

-0.37(89) 

ei5 

-10.0(25) 

-11.0(25) 

-11.0(25) 

-10.6(24) 

-10.5(24) 

-10.4(24) 

ei6 

+7.4(42) 

+7.0(42) 

+6.7(42) 

+4.6(42) 

+4.1(42) 

+3.5(41) 

cir 

+1.2(18) 

+1.2(18) 

+1.1(18) 

+1.3(18) 

+1.4(18) 

+1.4(18) 

ClS 

+8.4(19) 

+3.3(18) 

+3.6(17) 

-0.6(17) 

-1.5(17) 

-2.5(16) 




Table X. NNLOsim (A 

= 475) 



LEG 

TiX = 125 

TiX = 158 

Ta"" = 191 

Ta"" = 224 

Ta"" = 257 

Ta"" = 290 


-0.1513(20) 

-0.1507(20) 

-0.1503(20) 

-0.1493(20) 

-0.1489(20) 

-0.1483(20) 

Mpp) 

-0.1506(20) 

-0.1500(20) 

-0.1496(20) 

-0.1485(20) 

-0.1481(20) 

-0.1475(20) 

p^(nn) 

-0.1512(21) 

-0.1506(21) 

-0.1502(21) 

-0.1491(20) 

-0.1486(20) 

-0.1481(20) 

C^So 

+2.492(51) 

+2.464(50) 

+2.470(49) 

+2.508(48) 

+2.524(48) 

+2.541(47) 

Css, 

-0.1673(23) 

-0.1722(22) 

-0.1727(22) 

-0.1743(22) 

-0.1746(22) 

-0.1743(21) 

Css, 

+0.639(52) 

+0.638(50) 

+0.644(50) 

+0.690(48) 

+0.707(48) 

+0.723(47) 

Cei 

+0.297(24) 

+0.443(18) 

+0.455(17) 

+0.570(14) 

+0.599(13) 

+0.627(12) 

Cspq 

+1.036(36) 

+1.139(33) 

+1.130(33) 

+1.161(32) 

+1.162(31) 

+1.167(31) 

Cl Pi 

+0.062(42) 

+0.299(39) 

+0.318(38) 

+0.489(35) 

+0.540(35) 

+0.583(34) 

Csp, 

-0.857(34) 

-0.913(33) 

-0.925(33) 

-0.946(32) 

-0.954(31) 

-0.967(31) 

Csp^ 

-0.784(33) 

-0.700(32) 

-0.703(31) 

-0.663(31) 

-0.656(31) 

-0.649(30) 

CD 

+0.404(78) 

-0.003(64) 

-0.003(60) 

-0.317(54) 

-0.391(53) 

-0.471(51) 

Ce 

-0.175(22) 

-0.301(23) 

-0.304(21) 

-0.440(20) 

-0.475(20) 

-0.515(19) 

Cl 

-0.75(30) 

-0.49(30) 

-0.42(30) 

-0.12(30) 

-0.01(30) 

+0.11(30) 

C2 

+3.0(11) 

+3.2(11) 

+3.4(11) 

+4.1(11) 

+4.4(11) 

+4.7(10) 

C3 

-4.30(15) 

-3.78(14) 

-3.80(14) 

-3.58(14) 

-3.54(14) 

-3.51(13) 

C4 

+1.91(20) 

+2.85(15) 

+2.85(13) 

+3.68(11) 

+3.877(99) 

+4.092(90) 

di+d2 

+4.46(13) 

+4.81(11) 

+4.82(11) 

+5.20(10) 

+5.288(98) 

+5.391(95) 

ds 

-4.48(23) 

-4.58(23) 

-4.59(23) 

-4.75(22) 

-4.80(22) 

-4.85(22) 

d5 

+0.02(14) 

-0.11(14) 

-0.11(14) 

-0.22(14) 

-0.24(14) 

-0.27(14) 

di4 — di5 

-9.69(28) 

-9.90(27) 

-9.90(27) 

-10.15(27) 

-10.21(27) 

-10.28(27) 

ei4 

+0.93(91) 

+0.39(90) 

+0.38(90) 

-0.08(89) 

-0.20(89) 

-0.32(89) 

ei5 

-10.0(25) 

-11.2(25) 

-10.9(25) 

-10.4(24) 

-10.1(24) 

-9.7(24) 

ei6 

+6.9(42) 

+6.8(42) 

+6.0(42) 

+3.7(42) 

+2.7(41) 

+1.4(41) 

ei7 

+1.1(18) 

+1.1(18) 

+1.1(18) 

+1.3(18) 

+1.4(18) 

+1.5(18) 

ClS 

+8.3(19) 

+4.1(17) 

+4.1(17) 

+0.1(17) 

-0.8(16) 

-1.9(16) 



24 


Table XL NNLOsim (A = 500) 


LEG 

TiX = 125 

TiX = 158 

Ta""=i9i 

Ta"" = 224 

Ta"" = 257 

Ta"" = 290 


-0.1507(20) 

-0.1503(20) 

-0.1497(20) 

-0.1488(20) 

-0.1481(20) 

-0.1474(20) 

Mpp) 

-0.1500(20) 

-0.1496(20) 

-0.1490(20) 

-0.1480(20) 

-0.1473(20) 

-0.1465(20) 

P^(nn) 

-0.1508(21) 

-0.1503(21) 

-0.1497(21) 

-0.1486(20) 

-0.1479(20) 

-0.1471(20) 


+2.488(50) 

+2.451(49) 

+2.465(49) 

+2.502(48) 

+2.524(47) 

+2.548(47) 

C^s, 

-0.1605(24) 

-0.1668(23) 

-0.1673(23) 

-0.1691(22) 

-0.1693(22) 

-0.1687(21) 

C^s, 

+0.608(51) 

+0.601(50) 

+0.615(49) 

+0.660(48) 

+0.682(48) 

+0.705(47) 

Cei 

+0.287(22) 

+0.418(17) 

+0.436(16) 

+0.542(13) 

+0.572(12) 

+0.597(11) 

C^Po 

+1.096(35) 

+1.170(33) 

+1.161(32) 

+1.169(31) 

+1.164(31) 

+1.161(31) 

Cl Pi 

+0.071(42) 

+0.281(39) 

+0.294(37) 

+0.445(35) 

+0.487(34) 

+0.520(33) 

Cspi 

-0.813(34) 

-0.882(33) 

-0.900(32) 

-0.923(31) 

-0.936(31) 

-0.955(31) 

Csp^ 

-0.772(32) 

-0.694(32) 

-0.699(31) 

-0.663(31) 

-0.659(30) 

-0.658(30) 

CD 

+0.450(74) 

+0.102(62) 

+0.091(59) 

-0.189(54) 

-0.256(52) 

-0.325(51) 

Ce 

-0.297(15) 

-0.357(17) 

-0.362(17) 

-0.460(17) 

-0.490(17) 

-0.521(17) 

Cl 

-0.66(30) 

-0.45(30) 

-0.36(30) 

-0.07(30) 

+0.07(30) 

+0.22(30) 

C2 

+3.2(11) 

+3.3(11) 

+3.6(11) 

+4.3(11) 

+4.7(11) 

+5.1(10) 

C3 

-4.23(15) 

-3.75(14) 

-3.78(14) 

-3.58(14) 

-3.57(13) 

-3.56(13) 

C4 

+1.97(19) 

+2.73(14) 

+2.78(12) 

+3.527(99) 

+3.727(92) 

+3.933(85) 

di+d2 

+4.48(12) 

+4.74(11) 

+4.77(10) 

+5.115(98) 

+5.215(96) 

+5.320(94) 

ds 

-4.48(23) 

-4.54(23) 

-4.57(23) 

-4.72(22) 

-4.77(22) 

-4.83(22) 

d5 

+0.02(14) 

-0.10(14) 

-0.10(14) 

-0.19(14) 

-0.22(14) 

-0.24(14) 

di4 — di5 

-9.69(28) 

-9.84(27) 

-9.86(27) 

-10.09(27) 

-10.16(27) 

-10.23(27) 

ei4 

+0.88(90) 

+0.43(90) 

+0.40(90) 

-0.02(89) 

-0.14(89) 

-0.26(89) 

ei5 

-10.1(25) 

-11.4(25) 

-10.9(25) 

-10.4(24) 

-9.9(24) 

-9.3(24) 

ei6 

+6.3(42) 

+6.7(42) 

+5.5(42) 

+3.1(42) 

+1.7(41) 

-0.0(41) 

cir 

+1.1(18) 

+1.0(18) 

+1.1(18) 

+1.3(18) 

+1.4(18) 

+1.5(18) 

ClS 

+8.1(18) 

+4.7(17) 

+4.4(17) 

+0.9(16) 

-0.1(16) 

-1.2(16) 




Table XII. NNLOsim (A = 

525) 



LEG 

TiX = 125 

TiX = 158 

Ta"" = 191 

Ta"" = 224 

Ta"" = 257 

Ta"" = 290 


-0.1502(20) 

-0.1499(20) 

-0.1493(20) 

-0.1484(20) 

-0.1477(20) 

-0.1467(20) 

Mpp) 

-0.1495(20) 

-0.1492(20) 

-0.1485(20) 

-0.1476(20) 

-0.1468(20) 

-0.1459(20) 

p^(nn) 

-0.1504(21) 

-0.1500(21) 

-0.1493(21) 

-0.1483(20) 

-0.1475(20) 

-0.1466(20) 

C^So 

+2.484(50) 

+2.441(49) 

+2.460(49) 

+2.492(48) 

+2.517(47) 

+2.545(47) 

Css, 

-0.1540(26) 

-0.1616(24) 

-0.1625(23) 

-0.1646(22) 

-0.1648(22) 

-0.1639(22) 

Css, 

+0.581(51) 

+0.570(50) 

+0.590(49) 

+0.629(48) 

+0.655(47) 

+0.681(47) 

Cei 

+0.283(22) 

+0.399(17) 

+0.422(15) 

+0.519(12) 

+0.548(12) 

+0.571(11) 

Cspq 

+1.176(34) 

+1.219(33) 

+1.209(32) 

+1.193(31) 

+1.182(31) 

+1.171(30) 

Cl Pi 

+0.087(41) 

+0.273(38) 

+0.283(37) 

+0.416(34) 

+0.451(34) 

+0.475(33) 

Csp, 

-0.759(34) 

-0.844(33) 

-0.866(32) 

-0.890(31) 

-0.906(31) 

-0.930(31) 

Csp^ 

-0.760(32) 

-0.688(31) 

-0.694(31) 

-0.660(30) 

-0.659(30) 

-0.661(30) 

CD 

+0.508(72) 

+0.213(62) 

+0.193(59) 

-0.052(54) 

-0.111(53) 

-0.166(52) 

Ce 

-0.473(15) 

-0.462(14) 

-0.464(13) 

-0.521(14) 

-0.543(14) 

-0.566(14) 

Cl 

-0.59(30) 

-0.42(30) 

-0.31(30) 

-0.05(30) 

+0.10(30) 

+0.28(29) 

C2 

+3.3(11) 

+3.4(11) 

+3.7(11) 

+4.3(11) 

+4.8(10) 

+5.3(10) 

C3 

-4.17(15) 

-3.72(14) 

-3.75(14) 

-3.57(13) 

-3.57(13) 

-3.58(13) 

C4 

+2.06(17) 

+2.65(13) 

+2.73(11) 

+3.393(93) 

+3.589(88) 

+3.781(81) 

di+d2 

+4.50(12) 

+4.69(11) 

+4.74(10) 

+5.041(96) 

+5.143(95) 

+5.245(93) 

ds 

-4.49(23) 

-4.51(23) 

-4.55(23) 

-4.68(22) 

-4.74(22) 

-4.80(22) 

d5 

+0.01(14) 

-0.08(14) 

-0.09(14) 

-0.17(14) 

-0.19(14) 

-0.21(14) 

di4 — di5 

-9.70(28) 

-9.79(27) 

-9.83(27) 

-10.03(27) 

-10.10(27) 

-10.18(27) 

ei4 

+0.82(90) 

+0.46(90) 

+0.41(90) 

+0.03(89) 

-0.09(89) 

-0.20(89) 

ei5 

-10.2(25) 

-11.5(25) 

-10.9(25) 

-10.5(24) 

-9.9(24) 

-9.1(24) 

ei6 

+5.9(42) 

+6.5(42) 

+5.1(42) 

+3.0(41) 

+1.4(41) 

-0.7(41) 

ei7 

+1.1(18) 

+1.0(18) 

+1.1(18) 

+1.2(18) 

+1.4(18) 

+1.5(18) 

ClS 

+7.7(18) 

+5.1(17) 

+4.7(17) 

+1.5(16) 

+0.5(16) 

-0.4(16) 
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Table XIII. NNLOsim (A = 550) 


LEG 

TiX = 125 

TiX = 158 

Ta"" = 191 

Ta"" = 224 

Ta"" = 257 

Ta"" = 290 


-0.1498(20) 

-0.1495(20) 

-0.1489(20) 

-0.1482(20) 

-0.1474(20) 

-0.1464(20) 

Mpp) 

-0.1491(20) 

-0.1488(20) 

-0.1481(20) 

-0.1473(20) 

-0.1465(20) 

-0.1455(20) 

P^(nn) 

-0.1501(21) 

-0.1497(21) 

-0.1490(21) 

-0.1482(20) 

-0.1474(20) 

-0.1463(20) 


+2.480(50) 

+2.433(49) 

+2.454(48) 

+2.478(47) 

+2.504(47) 

+2.533(46) 

C^s, 

-0.1476(27) 

-0.1568(24) 

-0.1580(24) 

-0.1608(23) 

-0.1610(22) 

-0.1601(22) 

C^s, 

+0.557(51) 

+0.544(49) 

+0.567(49) 

+0.598(48) 

+0.626(47) 

+0.652(46) 

Cei 

+0.280(21) 

+0.385(16) 

+0.411(14) 

+0.500(12) 

+0.529(11) 

+0.550(11) 

C^Po 

+1.282(35) 

+1.292(33) 

+1.281(33) 

+1.237(31) 

+1.220(31) 

+1.203(30) 

Cl Pi 

+0.106(41) 

+0.274(38) 

+0.281(37) 

+0.400(34) 

+0.429(34) 

+0.447(33) 

Cspi 

-0.689(35) 

-0.796(33) 

-0.821(33) 

-0.844(31) 

-0.863(31) 

-0.890(31) 

Csp^ 

-0.748(32) 

-0.682(31) 

-0.687(31) 

-0.654(30) 

-0.654(30) 

-0.658(30) 

CD 

+0.587(71) 

+0.334(62) 

+0.307(60) 

+0.094(56) 

+0.043(55) 

+0.000(54) 

Ce 

-0.744(33) 

-0.646(19) 

-0.638(18) 

-0.645(15) 

-0.658(15) 

-0.673(14) 

Cl 

-0.54(30) 

-0.39(30) 

-0.28(30) 

-0.06(30) 

+0.09(30) 

+0.27(29) 

C2 

+3.4(11) 

+3.4(11) 

+3.8(11) 

+4.3(10) 

+4.7(10) 

+5.3(10) 

C3 

-4.10(15) 

-3.69(14) 

-3.72(14) 

-3.54(13) 

-3.55(13) 

-3.56(13) 

C4 

+2.15(17) 

+2.60(12) 

+2.70(11) 

+3.281(90) 

+3.467(84) 

+3.644(78) 

di+d2 

+4.54(12) 

+4.66(10) 

+4.72(10) 

+4.975(96) 

+5.073(94) 

+5.170(92) 

ds 

-4.49(23) 

-4.49(22) 

-4.53(22) 

-4.64(22) 

-4.70(22) 

-4.76(22) 

d5 

-0.01(14) 

-0.08(14) 

-0.08(14) 

-0.16(14) 

-0.18(14) 

-0.19(14) 

di4 — di5 

-9.72(27) 

-9.76(27) 

-9.81(27) 

-9.98(27) 

-10.05(27) 

-10.12(27) 

ei4 

+0.76(90) 

+0.47(90) 

+0.41(90) 

+0.08(89) 

-0.03(89) 

-0.14(89) 

ei5 

-10.3(25) 

-11.6(25) 

-11.0(24) 

-10.8(24) 

-10.1(24) 

-9.3(24) 

ei6 

+5.7(42) 

+6.4(42) 

+5.0(42) 

+3.4(41) 

+1.6(41) 

-0.5(41) 

cir 

+1.1(18) 

+1.0(18) 

+1.1(18) 

+1.2(18) 

+1.3(18) 

+1.4(18) 

ClS 

+7.3(18) 

+5.3(17) 

+4.8(17) 

+2.0(16) 

+1.1(16) 

+0.2(16) 




Table XIV. NNLOsim (A 

= 575) 



LEG 

TiX = 125 

TiX = 158 

Ta""=i9i 

Ta"" = 224 

Ta"" = 257 

Ta"" = 290 


-0.1495(20) 

-0.1492(20) 

-0.1486(20) 

-0.1481(20) 

-0.1474(20) 

-0.1464(20) 

Mpp) 

-0.1487(20) 

-0.1483(20) 

-0.1478(20) 

-0.1472(20) 

-0.1464(20) 

-0.1454(20) 

p^(nn) 

-0.1498(21) 

-0.1494(21) 

-0.1488(21) 

-0.1482(20) 

-0.1474(20) 

-0.1464(20) 

C^So 

+2.475(49) 

+2.426(49) 

+2.446(48) 

+2.461(47) 

+2.485(46) 

+2.512(46) 

Css, 

-0.1413(30) 

-0.1522(25) 

-0.1539(24) 

-0.1575(23) 

-0.1579(23) 

-0.1571(22) 

Css, 

+0.534(51) 

+0.522(49) 

+0.546(49) 

+0.568(47) 

+0.594(46) 

+0.619(46) 

Cei 

+0.280(21) 

+0.375(16) 

+0.404(14) 

+0.486(12) 

+0.513(11) 

+0.533(10) 

Cspq 

+1.426(37) 

+1.397(35) 

+1.386(34) 

+1.308(32) 

+1.284(31) 

+1.260(31) 

Cl Pi 

+0.127(41) 

+0.282(39) 

+0.286(37) 

+0.396(35) 

+0.421(34) 

+0.433(33) 

Csp, 

-0.598(37) 

-0.734(34) 

-0.762(33) 

-0.783(32) 

-0.804(31) 

-0.833(31) 

Csp^ 

-0.735(32) 

-0.676(31) 

-0.680(31) 

-0.645(30) 

-0.645(30) 

-0.649(30) 

CD 

+0.694(72) 

+0.472(64) 

+0.438(62) 

+0.251(59) 

+0.208(58) 

+0.176(57) 

Ce 

-1.232(78) 

-0.990(45) 

-0.955(40) 

-0.887(30) 

-0.887(28) 

-0.893(27) 

Cl 

-0.49(30) 

-0.37(30) 

-0.27(30) 

-0.10(30) 

+0.04(29) 

+0.21(29) 

C2 

+3.5(11) 

+3.5(11) 

+3.8(11) 

+4.1(10) 

+4.5(10) 

+5.1(10) 

C3 

-4.03(14) 

-3.66(14) 

-3.69(14) 

-3.50(13) 

-3.51(13) 

-3.52(13) 

C4 

+2.25(16) 

+2.57(12) 

+2.69(10) 

+3.188(87) 

+3.363(81) 

+3.523(76) 

di+d2 

+4.57(12) 

+4.63(10) 

+4.704(99) 

+4.914(95) 

+5.007(93) 

+5.096(92) 

ds 

-4.50(23) 

-4.47(22) 

-4.52(22) 

-4.60(22) 

-4.66(22) 

-4.72(22) 

d5 

-0.02(14) 

-0.07(14) 

-0.08(14) 

-0.15(14) 

-0.17(14) 

-0.18(14) 

di4 — di5 

-9.73(27) 

-9.74(27) 

-9.79(27) 

-9.93(27) 

-9.99(27) 

-10.06(27) 

ei4 

+0.70(90) 

+0.48(90) 

+0.41(89) 

+0.12(89) 

+0.02(89) 

-0.08(89) 

ei5 

-10.4(25) 

-11.8(25) 

-11.2(24) 

-11.2(24) 

-10.6(24) 

-9.8(24) 

ei6 

+5.5(42) 

+6.4(42) 

+5.1(42) 

+4.1(41) 

+2.5(40) 

+0.5(41) 

ei7 

+1.1(18) 

+1.0(18) 

+1.0(18) 

+1.1(18) 

+1.2(18) 

+1.3(18) 

ClS 

+6.9(17) 

+5.5(17) 

+4.9(16) 

+2.5(16) 

+1.7(16) 

+0.8(16) 
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Table XV. NNLOsim (A = 600) 


LEG 

TiX = 125 

TiX = 158 

Ta"" = 191 

Ta"" = 224 

Ta"" = 257 

Ta"" = 290 


-0.1491(20) 

-0.1487(20) 

-0.1483(20) 

-0.1481(20) 

-0.1474(20) 

-0.1466(20) 

Mpp) 

-0.1483(20) 

-0.1479(20) 

-0.1474(20) 

-0.1471(20) 

-0.1464(20) 

-0.1456(20) 

P^(nn) 

-0.1495(21) 

-0.1491(21) 

-0.1486(21) 

-0.1483(20) 

-0.1476(20) 

-0.1467(20) 


+2.469(49) 

+2.419(48) 

+2.437(48) 

+2.441(47) 

+2.461(46) 

+2.485(46) 

C^s, 

-0.1348(32) 

-0.1477(27) 

-0.1500(25) 

-0.1547(24) 

-0.1554(23) 

-0.1547(23) 

C^s, 

+0.512(50) 

+0.504(49) 

+0.527(48) 

+0.538(47) 

+0.563(47) 

+0.583(46) 

Cei 

+0.279(21) 

+0.368(15) 

+0.399(14) 

+0.476(12) 

+0.502(11) 

+0.520(10) 

C^Po 

+1.634(43) 

+1.556(38) 

+1.541(37) 

+1.416(33) 

+1.382(32) 

+1.351(31) 

Cl Pi 

+0.149(42) 

+0.296(40) 

+0.296(38) 

+0.400(35) 

+0.423(34) 

+0.432(33) 

Cspi 

-0.472(42) 

-0.653(35) 

-0.683(34) 

-0.703(32) 

-0.725(32) 

-0.756(31) 

Csp^ 

-0.723(32) 

-0.669(31) 

-0.671(31) 

-0.633(30) 

-0.632(30) 

-0.635(29) 

CD 

+0.835(74) 

+0.632(67) 

+0.592(65) 

+0.424(62) 

+0.388(62) 

+0.363(61) 

Ce 

-2.40(23) 

-1.76(13) 

-1.64(11) 

-1.409(77) 

-1.383(72) 

-1.372(68) 

Cl 

-0.46(30) 

-0.36(30) 

-0.28(30) 

-0.17(30) 

-0.04(29) 

+0.10(29) 

C2 

+3.5(11) 

+3.5(11) 

+3.7(11) 

+3.8(10) 

+4.2(10) 

+4.7(10) 

C3 

-3.97(14) 

-3.63(14) 

-3.65(14) 

-3.45(13) 

-3.45(13) 

-3.46(13) 

C4 

+2.33(15) 

+2.55(11) 

+2.69(10) 

+3.112(85) 

+3.274(80) 

+3.417(74) 

di+d2 

+4.60(11) 

+4.61(10) 

+4.692(99) 

+4.860(94) 

+4.944(93) 

+5.023(91) 

ds 

-4.51(23) 

-4.46(22) 

-4.51(22) 

-4.56(22) 

-4.61(22) 

-4.66(22) 

d5 

-0.03(14) 

-0.07(14) 

-0.09(14) 

-0.15(14) 

-0.16(14) 

-0.17(14) 

di4 — di5 

-9.75(27) 

-9.72(27) 

-9.78(27) 

-9.88(27) 

-9.94(27) 

-10.00(27) 

ei4 

+0.65(90) 

+0.48(90) 

+0.40(89) 

+0.16(90) 

+0.06(89) 

-0.03(89) 

ei5 

-10.5(25) 

-11.9(25) 

-11.5(24) 

-11.8(24) 

-11.3(24) 

-10.6(24) 

ei6 

+5.5(42) 

+6.4(42) 

+5.4(42) 

+5.2(41) 

+3.8(41) 

+2.1(40) 

cir 

+1.1(18) 

+0.9(18) 

+1.0(18) 

+1.1(18) 

+1.1(18) 

+1.2(18) 

ClS 

+6.5^7) 

+5.6(17) 

+4.9(16) 

+2.9(16) 

+2.1(16) 

+1.4(16) 



